module physics_msis

  use utils_constants
  use utils_spline
  use physics_constants

  private

  character(4) , public , dimension(2) :: chname , istime
  character(4) , public , dimension(3) :: isdate
  real(dp) , public , dimension(25) :: sav
  integer , public :: imr = 0
  real(dp) , public :: tlb , s , db04 , db16 , db28 , db32 , db40 , &
              db48 , db01 , za , t0 , z0 , xg0 , rl , dd , db14 , tr12

  public :: gtd7d

  public :: gtd7 , gts7 , ghp7
  public :: tselec , tretrv

  integer :: isw
  real(dp) , dimension(25) :: sw , swc
  real(dp) :: dm01 , dm04 , dm14 , dm16 , dm28 , dm32 , dm40
  real(dp) , dimension(10,8) :: pdm
  real(dp) , dimension(10) :: ptm
  real(dp) :: apd , apdf , c2tloc , c3tloc , ctloc , day , df , dfa ,   &
              s2tloc , s3tloc , stloc
  real(dp) , dimension(10) :: pavgm
  real(dp) , dimension(4) :: apt
  real(dp) , dimension(2) :: tgn1 , tgn2 , tgn3
  real(dp) , dimension(5) :: tn1 , tn3
  real(dp) , dimension(4) :: tn2
  integer :: iyr
  real(dp) , dimension(9,4) :: plg

  real(dp) , dimension(150) , target :: pt      ! pt1 , pt2 , pt3
  real(dp) , dimension(150,9) , target :: pd    ! pa1 , pa2 , pa3
                                                ! pb1 , pb2 , pb3
                                                ! pc1 , pc2 , pc3
                                                ! pd1 , pd2 , pd3
                                                ! pe1 , pe2 , pe3
                                                ! pf1 , pf2 , pf3
                                                ! pg1 , pg2 , pg3
                                                ! ph1 , ph2 , ph3
                                                ! pi1 , pi2 , pi3
  real(dp) , dimension(150) , target :: ps      ! pj1 , pj2 , pj3
  real(dp) , dimension(25,2) , target :: pdl    ! pk1
  real(dp) , dimension(100,4) , target :: ptl   ! pl1 , pl2
                                                ! pm1 , pm2
                                                ! pn1 , pn2
                                                ! po1 , po2
  real(dp) , dimension(100,10) , target :: pma  ! pp1 , pp2
                                                ! pq1 , pq2
                                                ! pr1 , pr2
                                                ! ps1 , ps2
                                                ! pu1 , pu2
                                                ! pv1 , pv2
                                                ! pw1 , pw2
                                                ! px1 , px2
                                                ! py1 , py2
                                                ! pz1 , pz2
  real(dp) , dimension(100) , target :: sam     ! paa1 , paa2

  real(dp) :: gsurf , re
  real(dp) :: gb , rout , tinf
  real(dp) , dimension(15) :: t
  real(dp) , parameter :: r100gas = rgasmol*d_100
  real(dp) , parameter :: nearzero = 0.000001_r8
  real(dp) , parameter :: dr = 1.72142E-2_r8
  integer , parameter :: mn1 = 5
  integer , parameter :: mn2 = 4
  integer , parameter :: mn3 = 5
  real(dp) , dimension(mn1) :: zn1
  real(dp) , dimension(mn2) :: zn2
  real(dp) , dimension(mn3) :: zn3
  integer , dimension(11) :: mt
  real(dp) , dimension(9) :: alpha
  real(dp) , dimension(8) :: altl

  ! MSISE-00 01-FEB-02
  data isdate/'01-F' , 'EB-0' , '2   '/
  data istime/'15:4' , '9:27'/
  data chname/'MSIS' , 'E-00'/
  !
  data zn1 /120.0_r8 , 110.0_r8 , 100.0_r8 , 90.0_r8 , 72.50_r8/
  data zn2 /72.5_r8 , 55.0_r8 , 45.0_r8 , 32.50_r8/
  data zn3 /32.5_r8 , 20.0_r8 , 15.0_r8 , 10.0_r8 , 0.0_r8/
  !
  data mt /48 , 0 , 4 , 16 , 28 , 32 , 40 , 1 , 49 , 14 , 17/
  data altl /200.0_r8 , 300.0_r8 , 160.0_r8 , 250.0_r8 , &
             240.0_r8 , 450.0_r8 , 320.0_r8 , 450.0_r8/
  data alpha / -0.380_r8 , 0.0_r8 , 0.0_r8 , 0.0_r8 , &
                0.170_r8 , 0.0_r8 , -0.380_r8 , 0.0_r8 , 0.0_r8/

  ! TEMPERATURE
  data pt(1:50) / &
            9.86573D-01 ,  1.62228D-02 ,  1.55270D-02 , -1.04323D-01 , &
           -3.75801D-03 , -1.18538D-03 , -1.24043D-01 ,  4.56820D-03 , &
            8.76018D-03 , -1.36235D-01 , -3.52427D-02 ,  8.84181D-03 , &
           -5.92127D-03 , -8.61650D+00 ,  0.00000D+00 ,  1.28492D-02 , &
            0.00000D+00 ,  1.30096D+02 ,  1.04567D-02 ,  1.65686D-03 , &
           -5.53887D-06 ,  2.97810D-03 ,  0.00000D+00 ,  5.13122D-03 , &
            8.66784D-02 ,  1.58727D-01 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 , -7.27026D-06 ,  0.00000D+00 ,  6.74494D+00 , &
            4.93933D-03 ,  2.21656D-03 ,  2.50802D-03 ,  0.00000D+00 , &
            0.00000D+00 , -2.08841D-02 , -1.79873D+00 ,  1.45103D-03 , &
            2.81769D-04 , -1.44703D-03 , -5.16394D-05 ,  8.47001D-02 , &
            1.70147D-01 ,  5.72562D-03 ,  5.07493D-05 ,  4.36148D-03 , &
            1.17863D-04 ,  4.74364D-03/
  data pt(51:100) / &
            6.61278D-03 ,  4.34292D-05 ,  1.44373D-03 ,  2.41470D-05 , &
            2.84426D-03 ,  8.56560D-04 ,  2.04028D-03 ,  0.00000D+00 , &
           -3.15994D+03 , -2.46423D-03 ,  1.13843D-03 ,  4.20512D-04 , &
            0.00000D+00 , -9.77214D+01 ,  6.77794D-03 ,  5.27499D-03 , &
            1.14936D-03 ,  0.00000D+00 , -6.61311D-03 , -1.84255D-02 , &
           -1.96259D-02 ,  2.98618D+04 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  6.44574D+02 ,  8.84668D-04 ,  5.05066D-04 , &
            0.00000D+00 ,  4.02881D+03 , -1.89503D-03 ,  0.00000D+00 , &
            0.00000D+00 ,  8.21407D-04 ,  2.06780D-03 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           -1.20410D-02 , -3.63963D-03 ,  9.92070D-05 , -1.15284D-04 , &
           -6.33059D-05 , -6.05545D-01 ,  8.34218D-03 , -9.13036D+01 , &
            3.71042D-04 ,  0.00000D+00/
  data pt(101:150) / &
            4.19000D-04 ,  2.70928D-03 ,  3.31507D-03 , -4.44508D-03 , &
           -4.96334D-03 , -1.60449D-03 ,  3.95119D-03 ,  2.48924D-03 , &
            5.09815D-04 ,  4.05302D-03 ,  2.24076D-03 ,  0.00000D+00 , &
            6.84256D-03 ,  4.66354D-04 ,  0.00000D+00 , -3.68328D-04 , &
            0.00000D+00 ,  0.00000D+00 , -1.46870D+02 ,  0.00000D+00 , &
            0.00000D+00 ,  1.09501D-03 ,  4.65156D-04 ,  5.62583D-04 , &
            3.21596D+00 ,  6.43168D-04 ,  3.14860D-03 ,  3.40738D-03 , &
            1.78481D-03 ,  9.62532D-04 ,  5.58171D-04 ,  3.43731D+00 , &
           -2.33195D-01 ,  5.10289D-04 ,  0.00000D+00 ,  0.00000D+00 , &
           -9.25347D+04 ,  0.00000D+00 , -1.99639D-03 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00/
  ! HE DENSITY
  data pd(1:50,1) / &
            1.09979D+00 , -4.88060D-02 , -1.97501D-01 ,                &
           -9.10280D-02 , -6.96558D-03 ,  2.42136D-02 ,  3.91333D-01 , &
           -7.20068D-03 , -3.22718D-02 ,  1.41508D+00 ,  1.68194D-01 , &
            1.85282D-02 ,  1.09384D-01 , -7.24282D+00 ,  0.00000D+00 , &
            2.96377D-01 , -4.97210D-02 ,  1.04114D+02 , -8.61108D-02 , &
           -7.29177D-04 ,  1.48998D-06 ,  1.08629D-03 ,  0.00000D+00 , &
            0.00000D+00 ,  8.31090D-02 ,  1.12818D-01 , -5.75005D-02 , &
           -1.29919D-02 , -1.78849D-02 , -2.86343D-06 ,  0.00000D+00 , &
           -1.51187D+02 , -6.65902D-03 ,  0.00000D+00 , -2.02069D-03 , &
            0.00000D+00 ,  0.00000D+00 ,  4.32264D-02 , -2.80444D+01 , &
           -3.26789D-03 ,  2.47461D-03 ,  0.00000D+00 ,  0.00000D+00 , &
            9.82100D-02 ,  1.22714D-01 , -3.96450D-02 ,  0.00000D+00 , &
           -2.76489D-03 ,  0.00000D+00 ,  1.87723D-03/
  data pd(51:100,1) / &
           -8.09813D-03 ,  4.34428D-05 , -7.70932D-03 ,                &
            0.00000D+00 , -2.28894D-03 , -5.69070D-03 , -5.22193D-03 , &
            6.00692D-03 , -7.80434D+03 , -3.48336D-03 , -6.38362D-03 , &
           -1.82190D-03 ,  0.00000D+00 , -7.58976D+01 , -2.17875D-02 , &
           -1.72524D-02 , -9.06287D-03 ,  0.00000D+00 ,  2.44725D-02 , &
            8.66040D-02 ,  1.05712D-01 ,  3.02543D+04 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 , -6.01364D+03 , -5.64668D-03 , &
           -2.54157D-03 ,  0.00000D+00 ,  3.15611D+02 , -5.69158D-03 , &
            0.00000D+00 ,  0.00000D+00 , -4.47216D-03 , -4.49523D-03 , &
            4.64428D-03 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  4.51236D-02 ,  2.46520D-02 ,  6.17794D-03 , &
            0.00000D+00 ,  0.00000D+00 , -3.62944D-01 , -4.80022D-02 , &
           -7.57230D+01 , -1.99656D-03 ,  0.00000D+00/
  data pd(101:150,1) / &
           -5.18780D-03 , -1.73990D-02 , -9.03485D-03 ,                &
            7.48465D-03 ,  1.53267D-02 ,  1.06296D-02 ,  1.18655D-02 , &
            2.55569D-03 ,  1.69020D-03 ,  3.51936D-02 , -1.81242D-02 , &
            0.00000D+00 , -1.00529D-01 , -5.10574D-03 ,  0.00000D+00 , &
            2.10228D-03 ,  0.00000D+00 ,  0.00000D+00 , -1.73255D+02 , &
            5.07833D-01 , -2.41408D-01 ,  8.75414D-03 ,  2.77527D-03 , &
           -8.90353D-05 , -5.25148D+00 , -5.83899D-03 , -2.09122D-02 , &
           -9.63530D-03 ,  9.77164D-03 ,  4.07051D-03 ,  2.53555D-04 , &
           -5.52875D+00 , -3.55993D-01 , -2.49231D-03 ,  0.00000D+00 , &
            0.00000D+00 ,  2.86026D+01 ,  0.00000D+00 ,  3.42722D-04 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00/
  ! O DENSITY
  data pd(1:50,2) / &
            1.02315D+00 , -1.59710D-01 , -1.06630D-01 ,                &
           -1.77074D-02 , -4.42726D-03 ,  3.44803D-02 ,  4.45613D-02 , &
           -3.33751D-02 , -5.73598D-02 ,  3.50360D-01 ,  6.33053D-02 , &
            2.16221D-02 ,  5.42577D-02 , -5.74193D+00 ,  0.00000D+00 , &
            1.90891D-01 , -1.39194D-02 ,  1.01102D+02 ,  8.16363D-02 , &
            1.33717D-04 ,  6.54403D-06 ,  3.10295D-03 ,  0.00000D+00 , &
            0.00000D+00 ,  5.38205D-02 ,  1.23910D-01 , -1.39831D-02 , &
            0.00000D+00 ,  0.00000D+00 , -3.95915D-06 ,  0.00000D+00 , &
           -7.14651D-01 , -5.01027D-03 ,  0.00000D+00 , -3.24756D-03 , &
            0.00000D+00 ,  0.00000D+00 ,  4.42173D-02 , -1.31598D+01 , &
           -3.15626D-03 ,  1.24574D-03 , -1.47626D-03 , -1.55461D-03 , &
            6.40682D-02 ,  1.34898D-01 , -2.42415D-02 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  6.13666D-04/
  data pd(51:100,2) / &
           -5.40373D-03 ,  2.61635D-05 , -3.33012D-03 ,                &
            0.00000D+00 , -3.08101D-03 , -2.42679D-03 , -3.36086D-03 , &
            0.00000D+00 , -1.18979D+03 , -5.04738D-02 , -2.61547D-03 , &
           -1.03132D-03 ,  1.91583D-04 , -8.38132D+01 , -1.40517D-02 , &
           -1.14167D-02 , -4.08012D-03 ,  1.73522D-04 , -1.39644D-02 , &
           -6.64128D-02 , -6.85152D-02 , -1.34414D+04 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  6.07916D+02 , -4.12220D-03 , &
           -2.20996D-03 ,  0.00000D+00 ,  1.70277D+03 , -4.63015D-03 , &
            0.00000D+00 ,  0.00000D+00 , -2.25360D-03 , -2.96204D-03 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  3.92786D-02 ,  1.31186D-02 , -1.78086D-03 , &
            0.00000D+00 ,  0.00000D+00 , -3.90083D-01 , -2.84741D-02 , &
           -7.78400D+01 , -1.02601D-03 ,  0.00000D+00/
  data pd(101:150,2) / &
           -7.26485D-04 , -5.42181D-03 , -5.59305D-03 ,                &
            1.22825D-02 ,  1.23868D-02 ,  6.68835D-03 , -1.03303D-02 , &
           -9.51903D-03 ,  2.70021D-04 , -2.57084D-02 , -1.32430D-02 , &
            0.00000D+00 , -3.81000D-02 , -3.16810D-03 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 , -9.05762D-04 , -2.14590D-03 , &
           -1.17824D-03 ,  3.66732D+00 , -3.79729D-04 , -6.13966D-03 , &
           -5.09082D-03 , -1.96332D-03 , -3.08280D-03 , -9.75222D-04 , &
            4.03315D+00 , -2.52710D-01 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00/
  ! N2 DENSITY
  data pd(1:50,3) / &
            1.16112D+00 ,  0.00000D+00 ,  0.00000D+00 ,  3.33725D-02 , &
            0.00000D+00 ,  3.48637D-02 , -5.44368D-03 ,  0.00000D+00 , &
           -6.73940D-02 ,  1.74754D-01 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  1.74712D+02 ,  0.00000D+00 ,  1.26733D-01 , &
            0.00000D+00 ,  1.03154D+02 ,  5.52075D-02 ,  0.00000D+00 , &
            0.00000D+00 ,  8.13525D-04 ,  0.00000D+00 ,  0.00000D+00 , &
            8.66784D-02 ,  1.58727D-01 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -2.50482D+01 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -2.48894D-03 , &
            6.16053D-04 , -5.79716D-04 ,  2.95482D-03 ,  8.47001D-02 , &
            1.70147D-01 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00/
  data pd(51:100,3) / &
           0.00000D+00 ,  2.47425D-05 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data pd(101:150,3) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  ! TLB
  data pd(1:50,4) / &
           9.44846D-01 ,  0.00000D+00 ,  0.00000D+00 , -3.08617D-02 , &
           0.00000D+00 , -2.44019D-02 ,  6.48607D-03 ,  0.00000D+00 , &
           3.08181D-02 ,  4.59392D-02 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  1.74712D+02 ,  0.00000D+00 ,  2.13260D-02 , &
           0.00000D+00 , -3.56958D+02 ,  0.00000D+00 ,  1.82278D-04 , &
           0.00000D+00 ,  3.07472D-04 ,  0.00000D+00 ,  0.00000D+00 , &
           8.66784D-02 ,  1.58727D-01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           3.83054D-03 ,  0.00000D+00 ,  0.00000D+00 , -1.93065D-03 , &
          -1.45090D-03 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -1.23493D-03 ,  1.36736D-03 ,  8.47001D-02 , &
           1.70147D-01 ,  3.71469D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data pd(51:100,4) / &
           5.10250D-03 ,  2.47425D-05 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  3.68756D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data pd(101:150,4) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  ! O2 DENSITY
  data pd(1:50,5) / &
           1.35580D+00 ,  1.44816D-01 ,  0.00000D+00 ,  6.07767D-02 , &
           0.00000D+00 ,  2.94777D-02 ,  7.46900D-02 ,  0.00000D+00 , &
          -9.23822D-02 ,  8.57342D-02 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  2.38636D+01 ,  0.00000D+00 ,  7.71653D-02 , &
           0.00000D+00 ,  8.18751D+01 ,  1.87736D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  1.49667D-02 ,  0.00000D+00 ,  0.00000D+00 , &
           8.66784D-02 ,  1.58727D-01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -3.67874D+02 , &
           5.48158D-03 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  8.47001D-02 , &
           1.70147D-01 ,  1.22631D-02 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data pd(51:100,5) / &
           8.17187D-03 ,  3.71617D-05 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -2.10826D-03 , -3.13640D-03 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
          -7.35742D-02 , -5.00266D-02 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  1.94965D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data pd(101:150,5) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  ! AR DENSITY
  data pd(1:50,6) / &
           1.04761D+00 ,  2.00165D-01 ,  2.37697D-01 ,  3.68552D-02 , &
           0.00000D+00 ,  3.57202D-02 , -2.14075D-01 ,  0.00000D+00 , &
          -1.08018D-01 , -3.73981D-01 ,  0.00000D+00 ,  3.10022D-02 , &
          -1.16305D-03 , -2.07596D+01 ,  0.00000D+00 ,  8.64502D-02 , &
           0.00000D+00 ,  9.74908D+01 ,  5.16707D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           8.66784D-02 ,  1.58727D-01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  3.46193D+02 , &
           1.34297D-02 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -3.48509D-03 , &
          -1.54689D-04 ,  0.00000D+00 ,  0.00000D+00 ,  8.47001D-02 , &
           1.70147D-01 ,  1.47753D-02 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data pd(51:100,6) / &
           1.89320D-02 ,  3.68181D-05 ,  1.32570D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  3.59719D-03 ,  7.44328D-03 , -1.00023D-03 , &
          -6.50528D+03 ,  0.00000D+00 ,  1.03485D-02 , -1.00983D-03 , &
          -4.06916D-03 , -6.60864D+01 , -1.71533D-02 ,  1.10605D-02 , &
           1.20300D-02 , -5.20034D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -2.62769D+03 ,  7.13755D-03 ,  4.17999D-03 , &
           0.00000D+00 ,  1.25910D+04 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -2.23595D-03 ,  4.60217D-03 ,  5.71794D-03 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
          -3.18353D-02 , -2.35526D-02 , -1.36189D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  2.03522D-02 , -6.67837D+01 , &
          -1.09724D-03 ,  0.00000D+00/
  data pd(101:150,6) / &
           -1.38821D-02 ,  1.60468D-02 ,  0.00000D+00 ,                &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  1.51574D-02 , -5.44470D-04 , &
            0.00000D+00 ,  7.28224D-02 ,  6.59413D-02 ,  0.00000D+00 , &
           -5.15692D-03 ,  0.00000D+00 ,  0.00000D+00 , -3.70367D+03 , &
            0.00000D+00 ,  0.00000D+00 ,  1.36131D-02 ,  5.38153D-03 , &
            0.00000D+00 ,  4.76285D+00 , -1.75677D-02 ,  2.26301D-02 , &
            0.00000D+00 ,  1.76631D-02 ,  4.77162D-03 ,  0.00000D+00 , &
            5.39354D+00 ,  0.00000D+00 , -7.51710D-03 ,  0.00000D+00 , &
            0.00000D+00 , -8.82736D+01 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00/
  ! H DENSITY
  data pd(1:50,7) / &
           1.26376D+00 , -2.14304D-01 , -1.49984D-01 ,  2.30404D-01 , &
           2.98237D-02 ,  2.68673D-02 ,  2.96228D-01 ,  2.21900D-02 , &
          -2.07655D-02 ,  4.52506D-01 ,  1.20105D-01 ,  3.24420D-02 , &
           4.24816D-02 , -9.14313D+00 ,  0.00000D+00 ,  2.47178D-02 , &
          -2.88229D-02 ,  8.12805D+01 ,  5.10380D-02 , -5.80611D-03 , &
           2.51236D-05 , -1.24083D-02 ,  0.00000D+00 ,  0.00000D+00 , &
           8.66784D-02 ,  1.58727D-01 , -3.48190D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  2.89885D-05 ,  0.00000D+00 ,  1.53595D+02 , &
          -1.68604D-02 ,  0.00000D+00 ,  1.01015D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  2.84552D-04 , &
          -1.22181D-03 ,  0.00000D+00 ,  0.00000D+00 ,  8.47001D-02 , &
           1.70147D-01 , -1.04927D-02 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -5.91313D-03/
  data pd(51:100,7) / &
           -2.30501D-02 ,  3.14758D-05 ,  0.00000D+00 ,                &
            0.00000D+00 ,  1.26956D-02 ,  8.35489D-03 ,  3.10513D-04 , &
            0.00000D+00 ,  3.42119D+03 , -2.45017D-03 , -4.27154D-04 , &
            5.45152D-04 ,  1.89896D-03 ,  2.89121D+01 , -6.49973D-03 , &
           -1.93855D-02 , -1.48492D-02 ,  0.00000D+00 , -5.10576D-02 , &
            7.87306D-02 ,  9.51981D-02 , -1.49422D+04 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  2.65503D+02 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  6.37110D-03 ,  3.24789D-04 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  6.14274D-02 ,  1.00376D-02 , -8.41083D-04 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -1.27099D-02 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00/
  data pd(101:150,7) / &
           -3.94077D-03 , -1.28601D-02 , -7.97616D-03 ,                &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 , -6.71465D-03 , -1.69799D-03 , &
            1.93772D-03 ,  3.81140D+00 , -7.79290D-03 , -1.82589D-02 , &
           -1.25860D-02 , -1.04311D-02 , -3.02465D-03 ,  2.43063D-03 , &
            3.63237D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00/
  ! N DENSITY
  data pd(1:50,8) / &
           7.09557D+01 , -3.26740D-01 ,  0.00000D+00 , -5.16829D-01 , &
          -1.71664D-03 ,  9.09310D-02 , -6.71500D-01 , -1.47771D-01 , &
          -9.27471D-02 , -2.30862D-01 , -1.56410D-01 ,  1.34455D-02 , &
          -1.19717D-01 ,  2.52151D+00 ,  0.00000D+00 , -2.41582D-01 , &
           5.92939D-02 ,  4.39756D+00 ,  9.15280D-02 ,  4.41292D-03 , &
           0.00000D+00 ,  8.66807D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           8.66784D-02 ,  1.58727D-01 ,  9.74701D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  6.70217D+01 , &
          -1.31660D-03 ,  0.00000D+00 , -1.65317D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  8.50247D-02 ,  2.77428D+01 ,  4.98658D-03 , &
           6.15115D-03 ,  9.50156D-03 , -2.12723D-02 ,  8.47001D-02 , &
           1.70147D-01 , -2.38645D-02 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  1.37380D-03/
  data pd(51:100,8) / &
           -8.41918D-03 ,  2.80145D-05 ,  7.12383D-03 ,                &
            0.00000D+00 , -1.66209D-02 ,  1.03533D-04 , -1.68898D-02 , &
            0.00000D+00 ,  3.64526D+03 ,  0.00000D+00 ,  6.54077D-03 , &
            3.69130D-04 ,  9.94419D-04 ,  8.42803D+01 , -1.16124D-02 , &
           -7.74414D-03 , -1.68844D-03 ,  1.42809D-03 , -1.92955D-03 , &
            1.17225D-01 , -2.41512D-02 ,  1.50521D+04 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  1.60261D+03 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 , -3.54403D-04 , -1.87270D-02 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  2.76439D-02 ,  6.43207D-03 , -3.54300D-02 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -2.80221D-02 , &
            8.11228D+01 , -6.75255D-04 ,  0.00000D+00/
  data pd(101:150,8) / &
           -1.05162D-02 , -3.48292D-03 , -6.97321D-03 ,                &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 , -1.45546D-03 , -1.31970D-02 , &
           -3.57751D-03 , -1.09021D+00 , -1.50181D-02 , -7.12841D-03 , &
           -6.64590D-03 , -3.52610D-03 , -1.87773D-02 , -2.22432D-03 , &
           -3.93895D-01 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00/
  ! HOT O DENSITY
  data pd(1:50,9) / &
           6.04050D-02 ,  1.57034D+00 ,  2.99387D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -1.51018D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -8.61650D+00 ,  1.26454D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  5.50878D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           8.66784D-02 ,  1.58727D-01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           6.23881D-02 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  8.47001D-02 , &
           1.70147D-01 , -9.45934D-02 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data pd(51:100,9) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data pd(101:150,9) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  ! S PARAM
  data ps(1:50) / &
           9.56827D-01 ,  6.20637D-02 ,  3.18433D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  3.94900D-02 ,  0.00000D+00 ,  0.00000D+00 , &
          -9.24882D-03 , -7.94023D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  1.74712D+02 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  2.74677D-03 ,  0.00000D+00 ,  1.54951D-02 , &
           8.66784D-02 ,  1.58727D-01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -6.99007D-04 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  1.24362D-02 , -5.28756D-03 ,  8.47001D-02 , &
           1.70147D-01 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data ps(51:100) / &
           0.00000D+00 ,  2.47425D-05 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data ps(101:150) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  ! TURBO
  data pdl/1.09930D+00 ,  3.90631D+00 ,  3.07165D+00 ,  9.86161D-01 , &
           1.63536D+01 ,  4.63830D+00 ,  1.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  1.28840D+00 ,  3.10302D-02 , &
           1.18339D-01 ,  1.00000D+00 ,  7.00000D-01 ,  1.15020D+00 , &
           3.44689D+00 ,  1.28840D+00 ,  1.00000D+00 ,  1.08738D+00 , &
           1.22947D+00 ,  1.10016D+00 ,  7.34129D-01 ,  1.15241D+00 , &
           2.22784D+00 ,  7.95046D-01 ,  4.01612D+00 ,  4.47749D+00 , &
           1.23435D+02 , -7.60535D-02 ,  1.68986D-06 ,  7.44294D-01 , &
           1.03604D+00 ,  1.72783D+02 ,  1.15020D+00 ,  3.44689D+00 , &
          -7.46230D-01 ,  9.49154D-01/

  ! LOWER BOUNDARY
  data ptm/1.04130D+03 ,  3.86000D+02 ,  1.95000D+02 ,  1.66728D+01 , &
           2.13000D+02 ,  1.20000D+02 ,  2.40000D+02 ,  1.87000D+02 , &
          -2.00000D+00 ,  0.00000D+00/
!
  data pdm/2.45600D+07 ,  6.71072D-06 ,  1.00000D+02 ,  0.00000D+00 , &
           1.10000D+02 ,  1.00000D+01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  8.59400D+10 ,  1.00000D+00 , &
           1.05000D+02 , -8.00000D+00 ,  1.10000D+02 ,  1.00000D+01 , &
           9.00000D+01 ,  2.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           2.81000D+11 ,  0.00000D+00 ,  1.05000D+02 ,  2.80000D+01 , &
           2.89500D+01 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  3.30000D+10 ,  2.68270D-01 , &
           1.05000D+02 ,  1.00000D+00 ,  1.10000D+02 ,  1.00000D+01 , &
           1.10000D+02 , -1.00000D+01 ,  0.00000D+00 ,  0.00000D+00 , &
           1.33000D+09 ,  1.19615D-02 ,  1.05000D+02 ,  0.00000D+00 , &
           1.10000D+02 ,  1.00000D+01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  1.76100D+05 ,  1.00000D+00 , &
           9.50000D+01 , -8.00000D+00 ,  1.10000D+02 ,  1.00000D+01 , &
           9.00000D+01 ,  2.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           1.00000D+07 ,  1.00000D+00 ,  1.05000D+02 , -8.00000D+00 , &
           1.10000D+02 ,  1.00000D+01 ,  9.00000D+01 ,  2.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  1.00000D+06 ,  1.00000D+00 , &
           1.05000D+02 , -8.00000D+00 ,  5.50000D+02 ,  7.60000D+01 , &
           9.00000D+01 ,  2.00000D+00 ,  0.00000D+00 ,  4.00000D+03/
  ! TN1(2)
  data ptl(1:50,1) / &
           1.00858D+00 ,  4.56011D-02 , -2.22972D-02 , -5.44388D-02 , &
           5.23136D-04 , -1.88849D-02 ,  5.23707D-02 , -9.43646D-03 , &
           6.31707D-03 , -7.80460D-02 , -4.88430D-02 ,  0.00000D+00 , &
           0.00000D+00 , -7.60250D+00 ,  0.00000D+00 , -1.44635D-02 , &
          -1.76843D-02 , -1.21517D+02 ,  2.85647D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  6.31792D-04 ,  0.00000D+00 ,  5.77197D-03 , &
           8.66784D-02 ,  1.58727D-01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -8.90272D+03 , &
           3.30611D-03 ,  3.02172D-03 ,  0.00000D+00 , -2.13673D-03 , &
          -3.20910D-04 ,  0.00000D+00 ,  0.00000D+00 ,  2.76034D-03 , &
           2.82487D-03 , -2.97592D-04 , -4.21534D-03 ,  8.47001D-02 , &
           1.70147D-01 ,  8.96456D-03 ,  0.00000D+00 , -1.08596D-02 , &
           0.00000D+00 ,  0.00000D+00/
  data ptl(51:100,1) / &
           5.57917D-03 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  9.65405D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  2.00000D+00/
  ! TN1(3)
  data ptl(1:50,2) / &
           9.39664D-01 ,  8.56514D-02 , -6.79989D-03 ,  2.65929D-02 , &
          -4.74283D-03 ,  1.21855D-02 , -2.14905D-02 ,  6.49651D-03 , &
          -2.05477D-02 , -4.24952D-02 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  1.19148D+01 ,  0.00000D+00 ,  1.18777D-02 , &
          -7.28230D-02 , -8.15965D+01 ,  1.73887D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -1.44691D-02 ,  2.80259D-04 , &
           8.66784D-02 ,  1.58727D-01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  2.16584D+02 , &
           3.18713D-03 ,  7.37479D-03 ,  0.00000D+00 , -2.55018D-03 , &
          -3.92806D-03 ,  0.00000D+00 ,  0.00000D+00 , -2.89757D-03 , &
          -1.33549D-03 ,  1.02661D-03 ,  3.53775D-04 ,  8.47001D-02 , &
           1.70147D-01 , -9.17497D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data ptl(51:100,2) / &
           3.56082D-03 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -1.00902D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  2.00000D+00/
  ! TN1(4)
  data ptl(1:50,3) / &
           9.85982D-01 , -4.55435D-02 ,  1.21106D-02 ,  2.04127D-02 , &
          -2.40836D-03 ,  1.11383D-02 , -4.51926D-02 ,  1.35074D-02 , &
          -6.54139D-03 ,  1.15275D-01 ,  1.28247D-01 ,  0.00000D+00 , &
           0.00000D+00 , -5.30705D+00 ,  0.00000D+00 , -3.79332D-02 , &
          -6.24741D-02 ,  7.71062D-01 ,  2.96315D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  6.81051D-03 , -4.34767D-03 , &
           8.66784D-02 ,  1.58727D-01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  1.07003D+01 , &
          -2.76907D-03 ,  4.32474D-04 ,  0.00000D+00 ,  1.31497D-03 , &
          -6.47517D-04 ,  0.00000D+00 , -2.20621D+01 , -1.10804D-03 , &
          -8.09338D-04 ,  4.18184D-04 ,  4.29650D-03 ,  8.47001D-02 , &
           1.70147D-01 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data ptl(51:100,3) / &
           -4.04337D-03 ,  0.00000D+00 ,  0.00000D+00 ,                &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -9.52550D-04 , &
            8.56253D-04 ,  4.33114D-04 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  1.21223D-03 ,  2.38694D-04 ,  9.15245D-04 , &
            1.28385D-03 ,  8.67668D-04 , -5.61425D-06 ,  1.04445D+00 , &
            3.41112D+01 ,  0.00000D+00 , -8.40704D-01 , -2.39639D+02 , &
            7.06668D-01 , -2.05873D+01 , -3.63696D-01 ,  2.39245D+01 , &
            0.00000D+00 , -1.06657D-03 , -7.67292D-04 ,  1.54534D-04 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  2.00000D+00/
  ! TN1(5) TN2(1)
  data ptl(1:50,4) / &
           1.00320D+00 ,  3.83501D-02 , -2.38983D-03 ,  2.83950D-03 , &
           4.20956D-03 ,  5.86619D-04 ,  2.19054D-02 , -1.00946D-02 , &
          -3.50259D-03 ,  4.17392D-02 , -8.44404D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  4.96949D+00 ,  0.00000D+00 , -7.06478D-03 , &
          -1.46494D-02 ,  3.13258D+01 , -1.86493D-03 ,  0.00000D+00 , &
          -1.67499D-02 ,  0.00000D+00 ,  0.00000D+00 ,  5.12686D-04 , &
           8.66784D-02 ,  1.58727D-01 , -4.64167D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  4.37353D-03 , -1.99069D+02 , &
           0.00000D+00 , -5.34884D-03 ,  0.00000D+00 ,  1.62458D-03 , &
           2.93016D-03 ,  2.67926D-03 ,  5.90449D+02 ,  0.00000D+00 , &
           0.00000D+00 , -1.17266D-03 , -3.58890D-04 ,  8.47001D-02 , &
           1.70147D-01 ,  0.00000D+00 ,  0.00000D+00 ,  1.38673D-02 , &
           0.00000D+00 ,  0.00000D+00/
  data ptl(51:100,4) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  1.60571D-03 ,  6.28078D-04 , &
           5.05469D-05 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
          -1.57829D-03 , -4.00855D-04 ,  5.04077D-05 , -1.39001D-03 , &
          -2.33406D-03 , -4.81197D-04 ,  1.46758D+00 ,  6.20332D+00 , &
           0.00000D+00 ,  3.66476D-01 , -6.19760D+01 ,  3.09198D-01 , &
          -1.98999D+01 ,  0.00000D+00 , -3.29933D+02 ,  0.00000D+00 , &
          -1.10080D-03 , -9.39310D-05 ,  1.39638D-04 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  2.00000D+00/
  ! TN2(2)
  data pma(1:50,1) / &
           9.81637D-01 , -1.41317D-03 ,  3.87323D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -3.58707D-02 , -8.63658D-03 ,  0.00000D+00 , &
           0.00000D+00 , -2.02226D+00 ,  0.00000D+00 , -8.69424D-03 , &
          -1.91397D-02 ,  8.76779D+01 ,  4.52188D-03 ,  0.00000D+00 , &
           2.23760D-02 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -7.07572D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -4.11210D-03 ,  3.50060D+01 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -8.36657D-03 ,  1.61347D+01 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -1.45130D-02 , &
           0.00000D+00 ,  0.00000D+00/
  data pma(51:100,1) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  1.24152D-03 ,  6.43365D-04 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           1.33255D-03 ,  2.42657D-03 ,  1.60666D-03 , -1.85728D-03 , &
          -1.46874D-03 , -4.79163D-06 ,  1.22464D+00 ,  3.53510D+01 , &
           0.00000D+00 ,  4.49223D-01 , -4.77466D+01 ,  4.70681D-01 , &
           8.41861D+00 , -2.88198D-01 ,  1.67854D+02 ,  0.00000D+00 , &
           7.11493D-04 ,  6.05601D-04 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  2.00000D+00/
  ! TN2(3)
  data pma(1:50,2) / &
           1.00422D+00 , -7.11212D-03 ,  5.24480D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -5.28914D-02 , -2.41301D-02 ,  0.00000D+00 , &
           0.00000D+00 , -2.12219D+01 , -1.03830D-02 , -3.28077D-03 , &
           1.65727D-02 ,  1.68564D+00 , -6.68154D-03 ,  0.00000D+00 , &
           1.45155D-02 ,  0.00000D+00 ,  8.42365D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -4.34645D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  2.16780D-02 ,  0.00000D+00 , -1.38459D+02 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  7.04573D-03 , -4.73204D+01 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  1.08767D-02 , &
           0.00000D+00 ,  0.00000D+00/
  data pma(51:100,2) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -8.08279D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  5.21769D-04 , -2.27387D-04 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           3.26769D-03 ,  3.16901D-03 ,  4.60316D-04 , -1.01431D-04 , &
           1.02131D-03 ,  9.96601D-04 ,  1.25707D+00 ,  2.50114D+01 , &
           0.00000D+00 ,  4.24472D-01 , -2.77655D+01 ,  3.44625D-01 , &
           2.75412D+01 ,  0.00000D+00 ,  7.94251D+02 ,  0.00000D+00 , &
           2.45835D-03 ,  1.38871D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  2.00000D+00/
  ! TN2(4) TN3(1)
  data pma(1:50,3) / &
           1.01890D+00 , -2.46603D-02 ,  1.00078D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -6.70977D-02 , -4.02286D-02 ,  0.00000D+00 , &
           0.00000D+00 , -2.29466D+01 , -7.47019D-03 ,  2.26580D-03 , &
           2.63931D-02 ,  3.72625D+01 , -6.39041D-03 ,  0.00000D+00 , &
           9.58383D-03 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -1.85291D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  1.39717D+02 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  9.19771D-03 , -3.69121D+02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -1.57067D-02 , &
           0.00000D+00 ,  0.00000D+00/
  data pma(51:100,3) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -7.07265D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -2.92953D-03 , -2.77739D-03 , &
          -4.40092D-04 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           2.47280D-03 ,  2.95035D-04 , -1.81246D-03 ,  2.81945D-03 , &
           4.27296D-03 ,  9.78863D-04 ,  1.40545D+00 , -6.19173D+00 , &
           0.00000D+00 ,  0.00000D+00 , -7.93632D+01 ,  4.44643D-01 , &
          -4.03085D+02 ,  0.00000D+00 ,  1.15603D+01 ,  0.00000D+00 , &
           2.25068D-03 ,  8.48557D-04 , -2.98493D-04 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  2.00000D+00/
  ! TN3(2)
  data pma(1:50,4) / &
           9.75801D-01 ,  3.80680D-02 , -3.05198D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  3.85575D-02 ,  5.04057D-02 ,  0.00000D+00 , &
           0.00000D+00 , -1.76046D+02 ,  1.44594D-02 , -1.48297D-03 , &
          -3.68560D-03 ,  3.02185D+01 , -3.23338D-03 ,  0.00000D+00 , &
           1.53569D-02 ,  0.00000D+00 , -1.15558D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  4.89620D-03 ,  0.00000D+00 , &
           0.00000D+00 , -1.00616D-02 , -8.21324D-03 , -1.57757D+02 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  6.63564D-03 ,  4.58410D+01 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -2.51280D-02 , &
           0.00000D+00 ,  0.00000D+00/
  data pma(51:100,4) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  9.91215D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -8.73148D-04 , -1.29648D-03 , &
          -7.32026D-05 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
          -4.68110D-03 , -4.66003D-03 , -1.31567D-03 , -7.39390D-04 , &
           6.32499D-04 , -4.65588D-04 , -1.29785D+00 , -1.57139D+02 , &
           0.00000D+00 ,  2.58350D-01 , -3.69453D+01 ,  4.10672D-01 , &
           9.78196D+00 , -1.52064D-01 , -3.85084D+03 ,  0.00000D+00 , &
          -8.52706D-04 , -1.40945D-03 , -7.26786D-04 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  2.00000D+00/
  ! TN3(3)
  data pma(1:50,5) / &
           9.60722D-01 ,  7.03757D-02 , -3.00266D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  2.22671D-02 ,  4.10423D-02 ,  0.00000D+00 , &
           0.00000D+00 , -1.63070D+02 ,  1.06073D-02 ,  5.40747D-04 , &
           7.79481D-03 ,  1.44908D+02 ,  1.51484D-04 ,  0.00000D+00 , &
           1.97547D-02 ,  0.00000D+00 , -1.41844D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  5.77884D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  9.74319D-03 ,  0.00000D+00 , -2.88015D+03 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -4.44902D-03 , -2.92760D+01 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  2.34419D-02 , &
           0.00000D+00 ,  0.00000D+00/
  data pma(51:100,5) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  5.36685D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -4.65325D-04 , -5.50628D-04 , &
           3.31465D-04 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
          -2.06179D-03 , -3.08575D-03 , -7.93589D-04 , -1.08629D-04 , &
           5.95511D-04 , -9.05050D-04 ,  1.18997D+00 ,  4.15924D+01 , &
           0.00000D+00 , -4.72064D-01 , -9.47150D+02 ,  3.98723D-01 , &
           1.98304D+01 ,  0.00000D+00 ,  3.73219D+03 ,  0.00000D+00 , &
          -1.50040D-03 , -1.14933D-03 , -1.56769D-04 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  2.00000D+00/
  ! TN3(4)
  data pma(1:50,6) / &
           1.03123D+00 , -7.05124D-02 ,  8.71615D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -3.82621D-02 , -9.80975D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  2.89286D+01 ,  9.57341D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  8.66153D+01 ,  7.91938D-04 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  4.68917D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  7.86638D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  9.90827D-03 ,  0.00000D+00 ,  6.55573D+01 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -4.00200D+01 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  7.07457D-03 , &
           0.00000D+00 ,  0.00000D+00/
  data pma(51:100,6) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  5.72268D-03 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -2.04970D-04 ,  1.21560D-03 , &
          -8.05579D-06 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
          -2.49941D-03 , -4.57256D-04 , -1.59311D-04 ,  2.96481D-04 , &
          -1.77318D-03 , -6.37918D-04 ,  1.02395D+00 ,  1.28172D+01 , &
           0.00000D+00 ,  1.49903D-01 , -2.63818D+01 ,  0.00000D+00 , &
           4.70628D+01 , -2.22139D-01 ,  4.82292D-02 ,  0.00000D+00 , &
          -8.67075D-04 , -5.86479D-04 ,  5.32462D-04 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  2.00000D+00/
  ! TN3(5) SURFACE TEMP TSL
  data pma(1:50,7) / &
           1.00828D+00 , -9.10404D-02 , -2.26549D-02 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -2.32420D-02 , -9.08925D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  3.36105D+01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -1.24957D+01 , -5.87939D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  2.79765D+01 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  2.01237D+03 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -1.75553D-02 , &
           0.00000D+00 ,  0.00000D+00/
  data pma(51:100,7) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  3.29699D-03 ,  1.26659D-03 , &
           2.68402D-04 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           1.17894D-03 ,  1.48746D-03 ,  1.06478D-04 ,  1.34743D-04 , &
          -2.20939D-03 , -6.23523D-04 ,  6.36539D-01 ,  1.13621D+01 , &
           0.00000D+00 , -3.93777D-01 ,  2.38687D+03 ,  0.00000D+00 , &
           6.61865D+02 , -1.21434D-01 ,  9.27608D+00 ,  0.00000D+00 , &
           1.68478D-04 ,  1.24892D-03 ,  1.71345D-03 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  2.00000D+00/
  ! TGN3(2) SURFACE GRAD TSLG
  data pma(1:50,8) / &
           1.57293D+00 , -6.78400D-01 ,  6.47500D-01 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 , -7.62974D-02 , -3.60423D-01 ,  0.00000D+00 , &
           0.00000D+00 ,  1.28358D+02 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  4.68038D+01 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 , -1.67898D-01 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  2.90994D+04 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  3.15706D+01 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           0.00000D+00 ,  0.00000D+00/
  data pma(51:100,8) / &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  &
           0.00000D+00 ,  2.00000D+00/
  ! TGN2(1) TGN1(2)
  data pma(1:50,9) / &
            8.60028D-01 ,  3.77052D-01 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 , -1.17570D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  7.77757D-03 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  1.01024D+02 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  6.54251D+02 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00/
  data pma(51:100,9) / &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
           -1.56959D-02 ,  1.91001D-02 ,  3.15971D-02 ,  1.00982D-02 , &
           -6.71565D-03 ,  2.57693D-03 ,  1.38692D+00 ,  2.82132D-01 , &
            0.00000D+00 ,  0.00000D+00 ,  3.81511D+02 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  2.00000D+00/
  ! TGN3(1) TGN2(2)
  data pma(1:50,10) / &
            1.06029D+00 , -5.25231D-02 ,  3.73034D-01 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  3.31072D-02 , -3.88409D-01 ,  0.00000D+00 , &
            0.00000D+00 , -1.65295D+02 , -2.13801D-01 , -4.38916D-02 , &
           -3.22716D-01 , -8.82393D+01 ,  1.18458D-01 ,  0.00000D+00 , &
           -4.35863D-01 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 , -1.19782D-01 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  2.62229D+01 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 , -5.37443D+01 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , -4.55788D-01 , &
            0.00000D+00 ,  0.00000D+00/
  data pma(51:100,10) / &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  3.84009D-02 ,  3.96733D-02 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            5.05494D-02 ,  7.39617D-02 ,  1.92200D-02 , -8.46151D-03 , &
           -1.34244D-02 ,  1.96338D-02 ,  1.50421D+00 ,  1.88368D+01 , &
            0.00000D+00 ,  0.00000D+00 , -5.13114D+01 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            5.11923D-02 ,  3.61225D-02 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 ,  0.00000D+00 , &
            0.00000D+00 ,  2.00000D+00/
  ! SEMIANNUAL MULT SAM
  data sam(1:50) / &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00 , 1.00000D+00 , 1.00000D+00 , &
            1.00000D+00 , 1.00000D+00/
  data sam(51:100) / &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00 , 0.00000D+00 , 0.00000D+00 , &
            0.00000D+00 , 0.00000D+00/
  ! MIDDLE ATMOSPHERE AVERAGES
  data pavgm /2.61000D+02 ,  2.64000D+02 ,  2.29000D+02 ,  2.17000D+02 , &
              2.17000D+02 ,  2.23000D+02 ,  2.86760D+02 , -2.93940D+00 , &
              2.50000D+00 ,  0.00000D+00/
  contains
!
!   Chemistry/dissociation correction for msis models
!
    subroutine meters(lm)
      implicit none
      logical , intent(in) :: lm
      if ( lm ) imr = 1
    end subroutine meters

    real(dp) function ccor(alt,r,h1,zh)
      implicit none
      ! alt - Altitude
      ! r   - Target ratio
      ! h1  - Transition scale length
      ! zh  - Altitude of 1/2 r
      real(dp) , intent(in) :: alt , h1 , r , zh
      real(dp) , parameter :: echeck = 70.0_r8
      real(dp) :: e , ex
      e = (alt-zh)/h1
      if ( e > echeck ) then
        ccor = d_0
      else if ( e < -echeck ) then
        ccor = r
      else
        ex = dexp(e)
        ccor = r/(d_1+ex)
      end if
      ccor = dexp(ccor)
    end function ccor
!
!   O&O2 Chemistry/dissociation correction for msis models
!
    real(dp) function ccor2(alt,r,h1,zh,h2)
      implicit none
      ! alt - Altitude
      ! r   - Target ratio
      ! h1  - Transition scale length
      ! zh  - Altitude of 1/2 r
      real(dp) , intent(in) :: alt , h1 , h2 , r , zh
      real(dp) :: e1 , e2 , ex1 , ex2
      real(dp) , parameter :: echeck = 70.0_r8
!
      e1 = (alt-zh)/h1
      e2 = (alt-zh)/h2
      if ( e1 > echeck .or. e2 > echeck ) then
        ccor2 = d_0
      else if ( e1 < -echeck .and. e2 < -echeck ) then
        ccor2 = r
      else
        ex1 = dexp(e1)
        ex2 = dexp(e2)
        ccor2 = r/(d_1+d_half*(ex1+ex2))
      end if
      ccor2 = dexp(ccor2)
    end function ccor2
!
! Calculate Temperature and Density Profiles for lower atmos.
!
    real(dp) function densm(alt,d0,xm,tz)
      implicit none
!
      real(dp) , intent(in) :: alt , d0 , xm
      real(dp) , intent(out) :: tz
!
      real(dp) :: expl , gamm , glb , t1 , t2 , x , y , yd1 , &
                 yd2 , yi , z , z1 , z2 , zg , zgdif
      integer :: k , mn
      real(dp) , dimension(10) :: xs , y2out , ys

      densm = d0
      if ( alt <= zn2(1) ) then
        ! STRATOSPHERE/MESOSPHERE TEMPERATURE
        z = dmax1(alt,zn2(mn2))
        mn = mn2
        z1 = zn2(1)
        z2 = zn2(mn)
        t1 = tn2(1)
        t2 = tn2(mn)
        zg = zeta(z,z1)
        zgdif = zeta(z2,z1)
        ! Set up spline nodes
        do k = 1 , mn
          xs(k) = zeta(zn2(k),z1)/zgdif
          ys(k) = d_1/tn2(k)
        end do
        yd1 = -tgn2(1)/(t1*t1)*zgdif
        yd2 = -tgn2(2)/(t2*t2)*zgdif*((re+z2)/(re+z1))**d_2
        ! Calculate spline coefficients
        call spline(xs(1:mn),ys(1:mn),yd1,yd2,y2out(1:mn))
        x = zg/zgdif
        call splint(xs(1:mn),ys(1:mn),y2out(1:mn),x,y)
        ! Temperature at altitude
        tz = d_1/y
        if ( dabs(xm) > nearzero ) then
          ! CALCULATE STRATOSPHERE/MESOSPHERE DENSITY
          glb = gsurf/(d_1+z1/re)**d_2
          gamm = xm*glb*zgdif/r100gas
          ! Integrate temperature profile
          call splini(xs(1:mn),ys(1:mn),y2out(1:mn),x,yi)
          expl = gamm*yi
          if ( expl > 50.0_r8 ) expl = 50.0_r8
          ! Density at altitude
          densm = densm*(t1/tz)*dexp(-expl)
        end if
        if ( alt <= zn3(1) ) then
          ! TROPOSPHERE/STRATOSPHERE TEMPERATURE
          z = alt
          mn = mn3
          z1 = zn3(1)
          z2 = zn3(mn)
          t1 = tn3(1)
          t2 = tn3(mn)
          zg = zeta(z,z1)
          zgdif = zeta(z2,z1)
          ! Set up spline nodes
          do k = 1 , mn
            xs(k) = zeta(zn3(k),z1)/zgdif
            ys(k) = d_1/tn3(k)
          end do
          yd1 = -tgn3(1)/(t1*t1)*zgdif
          yd2 = -tgn3(2)/(t2*t2)*zgdif*((re+z2)/(re+z1))**d_2
          ! Calculate spline coefficients
          call spline(xs(1:mn),ys(1:mn),yd1,yd2,y2out(1:mn))
          x = zg/zgdif
          call splint(xs(1:mn),ys(1:mn),y2out(1:mn),x,y)
          ! temperature at altitude
          tz = d_1/y
          if ( dabs(xm) > nearzero ) then
            ! CALCULATE TROPOSPHERIC/STRATOSPHERE DENSITY
            glb = gsurf/(d_1+z1/re)**d_2
            gamm = xm*glb*zgdif/r100gas
            ! Integrate temperature profile
            call splini(xs(1:mn),ys(1:mn),y2out(1:mn),x,yi)
            expl = gamm*yi
            if ( expl > 50.0_r8 ) expl = 50.0_r8
            ! Density at altitude
            densm = densm*(t1/tz)*dexp(-expl)
          end if
        end if
      end if
      if ( dabs(xm) < nearzero ) densm = tz
      contains
        real(dp) function zeta(zz,zl)
          implicit none
          real(dp) , intent(in) :: zz , zl
          zeta = (zz-zl)*(re+zl)/(re+zz)
        end function zeta
    end function densm
!
!   Calculate Temperature and Density Profiles for MSIS models
!   New lower thermo polynomial 10/30/89
!
    real(dp) function densu(alt,dlb,t1,t2,xm,xalph,tz,zlb,s2)
      implicit none
!
      real(dp) , intent(in) :: xalph , alt , t1 , t2 , dlb , s2 , xm , zlb
      real(dp) , intent(out) :: tz
!
      real(dp) :: densa , dta , expl , gamm , gammo , glb , &
                 tt1 , tt2 , ta , tt , x , y , yd1 , yd2 , yi , z ,   &
                 z1 , z2 , za , zg , zg2 , zgdif
      integer :: k , mn
      real(dp) , dimension(5) :: xs , y2out , ys
!
      densu = d_1
!     Joining altitude of Bates and spline
      za = zn1(1)
      z = dmax1(alt,za)
!     Geopotential altitude difference from ZLB
      zg2 = zeta(z,zlb)
!     Bates temperature
      tt = t1 - (t1-t2)*dexp(-s2*zg2)
      ta = tt
      tz = tt
      densu = tz
      if ( alt < za ) then
        ! CALCULATE TEMPERATURE BELOW ZA
        ! Temperature gradient at ZA from Bates profile
        dta = (t1-ta)*s2*((re+zlb)/(re+za))**d_2
        tgn1(1) = dta
        tn1(1) = ta
        z = dmax1(alt,zn1(mn1))
        mn = mn1
        z1 = zn1(1)
        z2 = zn1(mn)
        tt1 = tn1(1)
        tt2 = tn1(mn)
        ! Geopotental difference from Z1
        zg = zeta(z,z1)
        zgdif = zeta(z2,z1)
        ! Set up spline nodes
        do k = 1 , mn
          xs(k) = zeta(zn1(k),z1)/zgdif
          ys(k) = d_1/tn1(k)
        end do
        ! End node derivatives
        yd1 = -tgn1(1)/(tt1*tt1)*zgdif
        yd2 = -tgn1(2)/(tt2*tt2)*zgdif*((re+z2)/(re+z1))**d_2
        ! Calculate spline coefficients
        call spline(xs,ys,yd1,yd2,y2out)
        x = zg/zgdif
        call splint(xs,ys,y2out,x,y)
        ! temperature at altitude
        tz = d_1/y
        densu = tz
      end if
      if ( dabs(xm) > nearzero ) then
        ! CALCULATE DENSITY ABOVE ZA
        glb = gsurf/(d_1+zlb/re)**d_2
        gammo = xm*glb/(s2*r100gas*t1)
        expl = dexp(-s2*gammo*zg2)
        if ( expl > 50.0_r8 .or. tt <= d_0 ) expl = 50.0_r8
        ! Density at altitude
        densa = dlb*(t2/tt)**(d_1+xalph+gammo)*expl
        densu = densa
        if ( alt<za ) then
          ! CALCULATE DENSITY BELOW ZA
          glb = gsurf/(d_1+z1/re)**d_2
          gamm = xm*glb*zgdif/r100gas
          ! integrate spline temperatures
          call splini(xs,ys,y2out,x,yi)
          expl = gamm*yi
          if ( expl > 50.0_r8 .or. tz <= d_0 ) expl = 50.0_r8
          ! Density at altitude
          densu = densu*(tt1/tz)**(d_1+xalph)*exp(-expl)
        end if
      end if
      contains
        real(dp) function zeta(zz,zl)
          implicit none
          real(dp) , intent(in) :: zz , zl
          zeta = (zz-zl)*(re+zl)/(re+zz)
        end function zeta
    end function densu
!
!   Turbopause correction for msis models
!
    real(dp) function dnet(dd,dm,zhm,xmm,xm)
      implicit none
!
!     DD - diffusive density
!     DM - full mixed density
!     ZHM - transition scale length
!     XMM - full mixed molecular weight
!     XM  - species molecular weight
!     DNET - combined density
!
      real(dp) , intent(in) :: dm , xm , xmm , zhm
      real(dp) , intent(inout) :: dd
!
      real(dp) :: a , ylog
!
      a = zhm/(xmm-xm)
      if ( dm <= d_0 .or. dd <= d_0 ) then
        if ( dabs(dd) < nearzero .and. dabs(dm) < nearzero ) dd = d_1
        if ( dabs(dm) < nearzero ) then
          dnet = dd
        else if ( dabs(dd) < nearzero ) then
          dnet = dm
        end if
      else
        ylog = a*dlog(dm/dd)
        if ( ylog < -d_10 ) then
          dnet = dd
        else if ( ylog > d_10 ) then
          dnet = dm
        else
          dnet = dd*(d_1+dexp(ylog))**(d_1/a)
        end if
      end if
    end function dnet
!
!   SET SWITCHES
!   SW FOR MAIN TERMS, SWC FOR CROSS TERMS
!
!   TO TURN ON AND OFF PARTICULAR VARIATIONS CALL TSELEC(SV),
!   WHERE SV IS A 25 ELEMENT ARRAY CONTAINING 0. FOR OFF, 1.
!   FOR ON, OR 2. FOR MAIN EFFECTS OFF BUT CROSS TERMS ON
!
!   To get current values of SW: CALL TRETRV(SW)
!
    subroutine tselec(sv)
      implicit none
!
      real(dp) , dimension(25) , intent(in) :: sv
!
      integer :: i

      do i = 1 , 25
        sav(i) = sv(i)
        sw(i) = dmod(sv(i),d_2)
        if ( dabs(sv(i)-d_1) < nearzero .or. &
             dabs(sv(i)-d_2) < nearzero ) then
          swc(i) = d_1
        else
          swc(i) = d_0
        end if
      end do
      isw = 64999
      return
    end subroutine tselec
!
    subroutine tretrv(svv)
      real(dp) , dimension(25) , intent(out) :: svv
      do i = 1 , 25
        svv(i) = sav(i)
      end do
    end subroutine tretrv
!
!   CALCULATE LATITUDE VARIABLE GRAVITY (GV) AND EFFECTIVE RADIUS (REFF)
!
    subroutine glatf(lat)
      implicit none
      real(dp) , intent(in) :: lat
      real(dp)  :: c2
      c2 = dcos(d_2*deg2rad*lat)
      gsurf = egrav*d_100*(d_1-0.0026373_r8*c2)
      re = d_2*gsurf/(3.085462E-6_r8+2.27E-9_r8*c2)*1.E-5_r8
    end subroutine glatf
!
    real(dp) function glob7s(glong,p)
      implicit none
!
      real(dp) , intent(in) :: glong
      real(dp) , dimension(:) , intent(out) :: p
!
      real(dp) :: t82 , tt , t71 , t72 , t81
      integer :: i
      real(dp) , dimension(14) :: t
!
!     VERSION OF GLOBE FOR LOWER ATMOSPHERE 10/26/99
!
      real(dp) , parameter :: pset = d_2
      real(dp) , save :: dayl
      real(dp) , save :: p32 , p18 , p14 , p39
      real(dp) , save :: cd14 , cd18 , cd32 , cd39
      data dayl /-1.0_r8/
      data p32  /-1000.0_r8/
      data p18  /-1000.0_r8/
      data p14  /-1000.0_r8/
      data p39  /-1000.0_r8/

      t(:) = d_0

      if ( dabs(p(100)) < nearzero ) p(100) = pset

      if ( dabs(day-dayl) > nearzero .or. dabs(p32-p(32)) > nearzero ) then
        cd32 = dcos(dr*(day-p(32)))
      end if
      if ( dabs(day-dayl) > nearzero .or. dabs(p18-p(18)) > nearzero ) then
        cd18 = dcos(d_2*dr*(day-p(18)))
      end if
      if ( dabs(day-dayl) > nearzero .or. dabs(p14-p(14)) > nearzero ) then
        cd14 = dcos(dr*(day-p(14)))
      end if
      if ( dabs(day-dayl) > nearzero .or. dabs(p39-p(39)) > nearzero ) then
        cd39 = dcos(d_2*dr*(day-p(39)))
      end if

      dayl = day
      p32 = p(32)
      p18 = p(18)
      p14 = p(14)
      p39 = p(39)
!
      t(1) = p(22)*dfa
      ! TIME INDEPENDENT
      t(2) = p(2)*plg(3,1) + p(3)*plg(5,1) + p(23)*plg(7,1) + &
             p(27)*plg(2,1) + p(15)*plg(4,1) + p(60)*plg(6,1)
      ! SYMMETRICAL ANNUAL
      t(3) = (p(19)+p(48)*plg(3,1)+p(30)*plg(5,1))*cd32
      ! SYMMETRICAL SEMIANNUAL
      t(4) = (p(16)+p(17)*plg(3,1)+p(31)*plg(5,1))*cd18
      ! ASYMMETRICAL ANNUAL
      t(5) = (p(10)*plg(2,1)+p(11)*plg(4,1)+p(21)*plg(6,1))*cd14
      ! ASYMMETRICAL SEMIANNUAL
      t(6) = (p(38)*plg(2,1))*cd39
      ! DIURNAL
      if ( dabs(sw(7)) > d_0 ) then
        t71 = p(12)*plg(3,2)*cd14*swc(5)
        t72 = p(13)*plg(3,2)*cd14*swc(5)
        t(7) = ((p(4)*plg(2,2)+p(5)*plg(4,2)+t71)*ctloc + &
                (p(7)*plg(2,2)+p(8)*plg(4,2)+t72)*stloc)
      end if
      ! SEMIDIURNAL
      if ( dabs(sw(8)) > d_0 ) then
        t81 = (p(24)*plg(4,3)+p(36)*plg(6,3))*cd14*swc(5)
        t82 = (p(34)*plg(4,3)+p(37)*plg(6,3))*cd14*swc(5)
        t(8) = ((p(6)*plg(3,3)+p(42)*plg(5,3)+t81)*c2tloc + &
                (p(9)*plg(3,3)+p(43)*plg(5,3)+t82)*s2tloc)
      end if
      ! TERDIURNAL
      if ( dabs(sw(14)) > d_0 ) then
        t(14) = p(40)*plg(4,4)*s3tloc + p(41)*plg(4,4)*c3tloc
      end if
      ! MAGNETIC ACTIVITY
      if ( dabs(sw(9)) > d_0 ) then
        if ( sw(9) > d_0 ) t(9) = apdf*(p(33)+p(46)*plg(3,1)*swc(2))
        if ( sw(9) < d_0 ) t(9) = (p(51)*apt(1)+p(97)*plg(3,1)*apt(1)*swc(2))
      end if
      ! LONGITUDINAL
      if ( dabs(sw(10)) > d_0 .and. dabs(sw(11)) > d_0 .and. &
           glong > -1000.0_r8 ) then
        t(11) = (d_1+plg(2,1)*(p(81)*swc(5)*dcos(dr*(day-p(82))) +      &
                               p(86)*swc(6)*dcos(d_2*dr*(day-p(87)))) + &
                               p(84)*swc(3)*dcos(dr*(day-p(85))) +      &
                               p(88)*swc(4)*dcos(d_2*dr*(day-p(89)))) * &
                ((p(65)*plg(3,2)+p(66)*plg(5,2)+p(67)*plg(7,2)+         &
                  p(75)*plg(2,2)+p(76)*plg(4,2)+p(77)*plg(6,2)) *       &
                  dcos(deg2rad*glong) +                                   & 
                 (p(91)*plg(3,2)+p(92)*plg(5,2)+p(93)*plg(7,2)+         &
                  p(78)*plg(2,2)+p(79)*plg(4,2)+p(80)*plg(6,2)) *       &
                  dsin(deg2rad*glong))
      end if
      tt = d_0
      do i = 1 , 14
        tt = tt + abs(sw(i))*t(i)
      end do
      glob7s = tt
    end function glob7s
!
!     FIND ALTITUDE OF PRESSURE SURFACE (PRESS) FROM GTD7 INPUT:
!     IYD - YEAR AND DAY AS YYDDD
!     SEC - UT(SEC)
!     GLAT - GEODETIC LATITUDE(DEG)
!     GLONG - GEODETIC LONGITUDE(DEG)
!     STL - LOCAL APPARENT SOLAR TIME(HRS)
!     F107A - 3 MONTH AVERAGE OF F10.7 FLUX
!     F107 - DAILY F10.7 FLUX FOR PREVIOUS DAY
!     AP - MAGNETIC INDEX(DAILY) OR WHEN SW(9)=-1. :
!     - ARRAY CONTAINING:
!     (1) DAILY AP
!     (2) 3 HR AP INDEX FOR CURRENT TIME
!     (3) 3 HR AP INDEX FOR 3 HRS BEFORE CURRENT TIME
!     (4) 3 HR AP INDEX FOR 6 HRS BEFORE CURRENT TIME
!     (5) 3 HR AP INDEX FOR 9 HRS BEFORE CURRENT TIME
!     (6) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 12 TO 33 HRS PRIOR
!     TO CURRENT TIME
!     (7) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 36 TO 59 HRS PRIOR
!     TO CURRENT TIME
!     PRESS - PRESSURE LEVEL(MB)
!     OUTPUT:
!     ALT - ALTITUDE(KM)
!     D(1) - HE NUMBER DENSITY(CM-3)
!     D(2) - O NUMBER DENSITY(CM-3)
!     D(3) - N2 NUMBER DENSITY(CM-3)
!     D(4) - O2 NUMBER DENSITY(CM-3)
!     D(5) - AR NUMBER DENSITY(CM-3)
!     D(6) - TOTAL MASS DENSITY(GM/CM3)
!     D(7) - H NUMBER DENSITY(CM-3)
!     D(8) - N NUMBER DENSITY(CM-3)
!     D(9) - HOT O NUMBER DENSITY(CM-3)
!     T(1) - EXOSPHERIC TEMPERATURE
!     T(2) - TEMPERATURE AT ALT
!
    subroutine ghp7(iyd,sec,alt,glat,glong,stl,f107a,f107,ap,d,t,press)
      implicit none
!
      real(dp) , intent(in) :: f107 , f107a , glat , glong , press , sec , stl
      real(dp) , intent(out) :: alt
      integer , intent(in) :: iyd
      real(dp) , intent(in) , dimension(7) :: ap
      real(dp) , intent(inout) , dimension(9) :: d
      real(dp) , intent(inout) , dimension(2) :: t
!
      real(dp) :: ca , cd , cl , cl2 , diff , g , p , pl , sh , &
                 xm , xn , z , zi
      integer :: iday , l

      real(dp) , parameter :: bm = 1.3806E-19_r8
      real(dp) , parameter :: test = 0.00043_r8
      integer , parameter :: ltest = 12
!
      pl = dlog10(press)
!
!     Initial altitude estimate
      if ( pl >= -5.0_r8 ) then
        if ( pl > 2.5_r8 ) zi = 18.06_r8*(3.00_r8-pl)
        if ( pl > 0.75_r8 .and. pl <= 2.5_r8 ) zi = 14.98_r8*(3.08_r8-pl)
        if ( pl > -1.0_r8 .and. pl <= 0.75_r8 ) zi = 17.8_r8*(2.72_r8-pl)
        if ( pl > -2.0_r8 .and. pl <= -1.0_r8 ) zi = 14.28_r8*(3.64_r8-pl)
        if ( pl > -4.0_r8 .and. pl <= -2.0_r8 ) zi = 12.72_r8*(4.32_r8-pl)
        if ( pl <= -4.0_r8 ) zi = 25.3*(.11-pl)
        iday = mod(iyd,1000)
        cl = glat / d_090
        cl2 = cl*cl
        if ( iday < 182 ) cd = d_1 - iday/91.25_r8
        if ( iday >= 182 ) cd = iday / 91.25_r8 - 3.0_r8
        ca = d_0
        if ( pl > -1.11_r8 .and. pl <= -0.23_r8 ) ca = d_1
        if ( pl > -0.23_r8 ) ca = (2.79_r8-pl)/(2.79_r8+0.23_r8)
        if ( pl <= -1.11_r8 .and. pl > -3.0_r8 ) &
          ca = (-2.93_r8-pl)/(-2.93_r8+1.11_r8)
        z = zi - 4.87_r8*cl*cd*ca - 1.64_r8*cl2*ca + 0.31_r8*ca*cl
      end if
      if ( pl < -5.0_r8 ) z = 22.0_r8*(pl+4.0_r8)**d_2 + 110.0_r8
      ! ITERATION LOOP
      l = 0
      do
        l = l + 1
        call gtd7(iyd,sec,z,glat,glong,stl,f107a,f107,ap,48,d,t)
        xn = d(1) + d(2) + d(3) + d(4) + d(5) + d(7) + d(8)
        p = bm*xn*t(2)
        if ( imr == 1 ) p = p*1.E-6_r8
        diff = pl - dlog10(p)
        if ( abs(diff) < test .or. l == ltest ) then
          ! Non converging , should not happen
          alt = z
          exit
        else
          xm = d(6)/xn/1.66E-24_r8
          if ( imr == 1 ) xm = xm*1.E3_r8
          g = gsurf / (1.0_r8+z/re)**d_2
          sh = r100gas*t(2)/(xm*g)
          ! New altitude estimate using scale height
          if ( l < 6 ) then
            z = z - sh*diff*2.302_r8
          else
            z = z - sh*diff
          end if
          cycle
        end if
      end do
    end subroutine ghp7
!
    real(dp) function globe7(yrd,sec,lat,long,tloc,f107a,f107,ap,p)
      implicit none
!
      real(dp) , intent(in) :: f107 , f107a , lat , long , sec , tloc , yrd
      real(dp) , dimension(:) , intent(in) :: ap
      real(dp) , intent(out) , dimension(:) :: p
!
      real(dp) :: c , c2 , c4 , exp1 , f1 , f2 , p44 , p45 , s , s2 ,    &
                  t71 , t72 , t81 , t82 , tix
      integer :: i
      real(dp) , save , dimension(25) :: sv = d_1
!
!     CALCULATE G(L) FUNCTION
!
!     Upper Thermosphere Parameters
      integer , parameter :: nsw = 14
      real(dp) , parameter :: hr = 0.2618_r8
      real(dp) , parameter :: sr = 7.2722E-5_r8
      real(dp) , save :: xl  = 1000.0_r8
      real(dp) , save :: tll = 1000.0_r8
      real(dp) , save :: sw9 = 1.0_r8
      real(dp) , save :: dayl = -1.0_r8
      real(dp) , save :: p14 = -1000.0_r8
      real(dp) , save :: p18 = -1000.0_r8
      real(dp) , save :: p32 = -1000.0_r8
      real(dp) , save :: p39 = -1000.0_r8
      real(dp) , save :: cd14 , cd18 , cd32 , cd39

      if ( isw /= 64999 ) call tselec(sv)
!
      t(:) = d_0
!
      if ( sw(9) > d_0 ) sw9 =  d_1
      if ( sw(9) < d_0 ) sw9 = -d_1

      iyr = idint(yrd*d_r1000)
      day = yrd - iyr*d_1000
!     Eq. A22 (remainder of code)
      if ( dabs(xl-lat) > nearzero ) then
!       CALCULATE LEGENDRE POLYNOMIALS
        c = dsin(lat*deg2rad)
        s = dcos(lat*deg2rad)
        c2 = c*c
        c4 = c2*c2
        s2 = s*s
        plg(2,1) = c
        plg(3,1) = 0.5_r8*(d_3*c2-d_1)
        plg(4,1) = 0.5_r8*(d_5*c*c2-d_3*c)
        plg(5,1) = (35.0_r8*c4-30.0_r8*c2+d_3)/d_8
        plg(6,1) = (63.0_r8*c2*c2*c-70.0_r8*c2*c+15.0_r8*c)/d_8
        plg(7,1) = (11.0_r8*c*plg(6,1)-d_5*plg(5,1))/d_6
!       PLG(8,1) = (13.0_r8*C*PLG(7,1) - d_6*PLG(6,1))/7.0_r8
        plg(2,2) = s
        plg(3,2) = d_3*c*s
        plg(4,2) = 1.5_r8*(d_5*c2-d_1)*s
        plg(5,2) = 2.5_r8*(d_7*c2*c-d_3*c)*s
        plg(6,2) = 1.875_r8*(21.0_r8*c4-14.0_r8*c2+d_1)*s
        plg(7,2) = (11.0_r8*c*plg(6,2)-d_6*plg(5,2))/d_5
!       PLG(8,2) = (13.0_r8*C*PLG(7,2)-d_7*PLG(6,2))/d_6
!       PLG(9,2) = (15.0_r8*C*PLG(8,2)-d_8*PLG(7,2))/d_7
        plg(3,3) = d_3*s2
        plg(4,3) = 15.0_r8*s2*c
        plg(5,3) = 7.5_r8*(d_7*c2-d_1)*s2
        plg(6,3) = d_3*c*plg(5,3) - d_2*plg(4,3)
        plg(7,3) = (11.0_r8*c*plg(6,3)-d_7*plg(5,3))/d_4
        plg(8,3) = (13.0_r8*c*plg(7,3)-d_8*plg(6,3))/d_5
        plg(4,4) = 15.0_r8*s2*s
        plg(5,4) = 105.0_r8*s2*s*c
        plg(6,4) = (d_9*c*plg(5,4)-d_7*plg(4,4))/d_2
        plg(7,4) = (11.0_r8*c*plg(6,4)-d_8*plg(5,4))/d_3
        xl = lat
      end if
      if ( dabs(tll-tloc) > nearzero ) then
        if ( dabs(sw(7)) > d_0 .or. dabs(sw(8)) > d_0 .or. &
             dabs(sw(14)) > d_0 ) then
          stloc = dsin(hr*tloc)
          ctloc = dcos(hr*tloc)
          s2tloc = dsin(d_2*hr*tloc)
          c2tloc = dcos(d_2*hr*tloc)
          s3tloc = dsin(d_3*hr*tloc)
          c3tloc = dcos(d_3*hr*tloc)
          tll = tloc
        end if
      end if
      if ( dabs(day-dayl) > nearzero .or. dabs(p(14)-p14) > nearzero ) &
        cd14 = dcos(dr*(day-p(14)))
      if ( dabs(day-dayl) > nearzero .or. dabs(p(18)-p18) > nearzero ) &
        cd18 = dcos(d_2*dr*(day-p(18)))
      if ( dabs(day-dayl) > nearzero .or. dabs(p(32)-p32) > nearzero ) &
        cd32 = dcos(dr*(day-p(32)))
      if ( dabs(day-dayl) > nearzero .or. dabs(p(39)-p39) > nearzero ) &
        cd39 = dcos(d_2*dr*(day-p(39)))
      dayl = day
      p14 = p(14)
      p18 = p(18)
      p32 = p(32)
      p39 = p(39)
      ! F10.7 EFFECT
      df = f107 - f107a
      dfa = f107a - 150.0_r8
      t(1) = p(20)*df*(d_1+p(60)*dfa) + p(21)*df*df + p(22)*dfa+p(30)*dfa*dfa
      f1 = d_1 + (p(48)*dfa+p(20)*df+p(21)*df*df)*swc(1)
      f2 = d_1 + (p(50)*dfa+p(20)*df+p(21)*df*df)*swc(1)
      ! TIME INDEPENDENT
      t(2) = (p(2)*plg(3,1)+p(3)*plg(5,1)+p(23)*plg(7,1)) + &
             (p(15)*plg(3,1))*dfa*swc(1) + p(27)*plg(2,1)
      ! SYMMETRICAL ANNUAL
      t(3) = (p(19))*cd32
      ! SYMMETRICAL SEMIANNUAL
      t(4) = (p(16)+p(17)*plg(3,1))*cd18
      ! ASYMMETRICAL ANNUAL
      t(5) = f1*(p(10)*plg(2,1)+p(11)*plg(4,1))*cd14
      ! ASYMMETRICAL SEMIANNUAL
      t(6) = p(38)*plg(2,1)*cd39
      ! DIURNAL
      if ( dabs(sw(7)) > d_0 ) then
        t71 = (p(12)*plg(3,2))*cd14*swc(5)
        t72 = (p(13)*plg(3,2))*cd14*swc(5)
        t(7) = f2*((p(4)*plg(2,2)+p(5)*plg(4,2)+p(28)*plg(6,2)+t71)*ctloc + &
                   (p(7)*plg(2,2)+p(8)*plg(4,2)+p(29)*plg(6,2)+t72)*stloc)
      end if
      ! SEMIDIURNAL
      if ( dabs(sw(8)) > d_0 ) then
        t81 = (p(24)*plg(4,3)+p(36)*plg(6,3))*cd14*swc(5)
        t82 = (p(34)*plg(4,3)+p(37)*plg(6,3))*cd14*swc(5)
        t(8) = f2*((p(6)*plg(3,3)+p(42)*plg(5,3)+t81)*c2tloc + &
                   (p(9)*plg(3,3)+p(43)*plg(5,3)+t82)*s2tloc)
      end if
      ! TERDIURNAL
      if ( dabs(sw(14)) > d_0 ) then
        t(14) = f2*((p(40)*plg(4,4) + &
              (p(94)*plg(5,4)+p(47)*plg(7,4))*cd14*swc(5))*s3tloc + &
              (p(41)*plg(4,4)+(p(95)*plg(5,4) + &
               p(49)*plg(7,4))*cd14*swc(5))*c3tloc)
      end if
      ! MAGNETIC ACTIVITY BASED ON DAILY AP
      if ( dabs(sw9+d_1) > nearzero ) then
        apd = (ap(1)-d_4)
        p44 = p(44)
        p45 = p(45)
        if ( p44 < d_0 ) p44 = 1.E-5_r8
        apdf = apd + (p45-d_1)*(apd+(dexp(-p44*apd)-d_1)/p44)
        if ( dabs(sw(9)) > nearzero ) then
          t(9) = apdf*(p(33)+p(46)*plg(3,1)+p(35)*plg(5,1) +           &
                      (p(101)*plg(2,1)+p(102)*plg(4,1) +               &
                       p(103)*plg(6,1))*cd14*swc(5)+(p(122)*plg(2,2) + &
                       p(123)*plg(4,2)+p(124)*plg(6,2))*swc(7)*        &
                       dcos(hr*(tloc-p(125))))
        end if
      else if ( dabs(p(52)) > nearzero ) then
        exp1 = dexp(-10800.0_r8*dabs(p(52))/(d_1+p(139)*(d_045-dabs(lat))))
        if ( exp1 > 0.99999_r8 ) exp1 = 0.99999_r8
        if ( p(25) < 1.E-4_r8 ) p(25) = 1.E-4_r8
        apt(1) = sg0(exp1)
!       APT(2)=SG2(EXP1)
!       APT(3)=SG0(EXP2)
!       APT(4)=SG2(EXP2)
        if ( dabs(sw(9)) > d_0 ) then
          t(9) = apt(1)*(p(51)+p(97)*plg(3,1)+p(55)*plg(5,1) +   &
                        (p(126)*plg(2,1)+p(127)*plg(4,1) +       &
                         p(128)*plg(6,1))*cd14*swc(5) +          &
                        (p(129)*plg(2,2)+p(130)*plg(4,2) +       &
                         p(131)*plg(6,2))*swc(7)*dcos(hr*(tloc-p(132))))
        end if
      end if
      if ( dabs(sw(10)) > d_0 .and. long > -1000.0_r8 ) then
        ! LONGITUDINAL
        if ( dabs(sw(11)) > d_0 ) then
          t(11) = (d_1+p(81)*dfa*swc(1)) *                               &
                    ((p(65)*plg(3,2)+p(66)*plg(5,2)+p(67)*plg(7,2) +     &
                      p(104)*plg(2,2)+p(105)*plg(4,2)+p(106)*plg(6,2) +  &
              swc(5)*(p(110)*plg(2,2)+p(111)*plg(4,2)+p(112)*plg(6,2)) * &
                  cd14)*dcos(deg2rad*long) +                             &
                     (p(91)*plg(3,2)+p(92)*plg(5,2)+p(93)*plg(7,2) +     &
                      p(107)*plg(2,2)+p(108)*plg(4,2)+p(109)*plg(6,2) +  &
              swc(5)*(p(113)*plg(2,2)+p(114)*plg(4,2)+p(115)*plg(6,2)) * &
                  cd14)*dsin(deg2rad*long))
        end if
        ! UT AND MIXED UT,LONGITUDE
        if ( dabs(sw(12)) > nearzero ) then
          t(12) = (d_1+p(96)*plg(2,1))*(d_1+p(82)*dfa*swc(1)) *   &
                  (d_1+p(120)*plg(2,1)*swc(5)*cd14)  *            &
                ((p(69)*plg(2,1)+p(70)*plg(4,1)+p(71)*plg(6,1)) * &
                dcos(sr*(sec-p(72))))
          t(12) = t(12) + swc(11) *                              &
                (p(77)*plg(4,3)+p(78)*plg(6,3)+p(79)*plg(8,3)) * &
                dcos(sr*(sec-p(80))+d_2*deg2rad*long)*(d_1+p(138)*dfa*swc(1))
        end if
        ! UT,LONGITUDE MAGNETIC ACTIVITY
        if ( dabs(sw(13)) > nearzero ) then
          if ( dabs(sw9+d_1) > nearzero ) then
            t(13) = apdf*swc(11)*(d_1+p(121)*plg(2,1))                 &
                   *((p(61)*plg(3,2)+p(62)*plg(5,2)+p(63)*plg(7,2))    &
                   *dcos(deg2rad*(long-p(64)))) + apdf*swc(11)*swc(5)  &
                   *(p(116)*plg(2,2)+p(117)*plg(4,2)+p(118)*plg(6,2))  &
                   *cd14*dcos(deg2rad*(long-p(119))) + apdf*swc(12)    &
                   *(p(84)*plg(2,1)+p(85)*plg(4,1)+p(86)*plg(6,1))     &
                   *dcos(sr*(sec-p(76)))
          else if ( p(52)/=0 ) then
            t(13) = apt(1)*swc(11)*(d_1+p(133)*plg(2,1))                &
                   *((p(53)*plg(3,2)+p(99)*plg(5,2)+p(68)*plg(7,2))     &
                   *dcos(deg2rad*(long-p(98)))) + apt(1)*swc(11)*swc(5) &
                   *(p(134)*plg(2,2)+p(135)*plg(4,2)+p(136)*plg(6,2))   &
                   *cd14*dcos(deg2rad*(long-p(137))) + apt(1)*swc(12)   &
                   *(p(56)*plg(2,1)+p(57)*plg(4,1)+p(58)*plg(6,1))      &
                   *dcos(sr*(sec-p(59)))
          end if
        end if
      end if

      ! PARMS NOT USED: 83, 90,100,140-150
      tix = p(31)
      do i = 1 , nsw
        tix = tix + abs(sw(i))*t(i)
      end do
      globe7 = tix
      contains

!     3hr Magnetic activity functions
!     Eq. A24d
      real(dp) function g0(a)
        implicit none
        real(dp) , intent(in) :: a
        g0 = (a-d_4+(p(26)-d_1)*(a-d_4 + &
             (dexp(-dabs(p(25))*(a-d_4))-d_1)/dabs(p(25))))
      end function g0
!     Eq. A24c
      real(dp) function sumex(ex)
        implicit none
        real(dp) , intent(in) :: ex
        sumex = d_1 + (d_1-ex**19.0_r8)/(d_1-ex)*ex**(d_half)
      end function sumex
!     Eq. A24a
      real(dp) function sg0(ex)
        implicit none
        real(dp) , intent(in) :: ex
        sg0 = (g0(ap(2))+(g0(ap(3))*ex+g0(ap(4))*ex*ex + &
               g0(ap(5))*ex**d_3+(g0(ap(6))*ex**d_4 +    &
               g0(ap(7))*ex**12.0_r8)*(d_1-ex**d_8)/(d_1-ex)))/sumex(ex)
      end function sg0
    end function globe7
!
    subroutine gtd7d(iyd,sec,alt,glat,glong,stl,f107a,f107,ap,mass,d,t)
      implicit none
!
      real(dp) , intent(in) :: alt , f107 , f107a , glat , glong , sec , stl
      integer :: iyd , mass
      real(dp) , intent(in) , dimension(7) :: ap
      real(dp) , intent(inout) , dimension(2) :: t
      real(dp) , intent(inout) , dimension(9) :: d
!
!     NRLMSISE-00
!     -----------
!     This subroutine provides Effective Total Mass Density for
!     output D(6) which includes contributions from "anomalous
!     oxygen" which can affect satellite drag above 500 km.  This
!     subroutine is part of the distribution package for the
!     Neutral Atmosphere Empirical Model from the surface to lower
!     exosphere.  See subroutine GTD7 for more extensive comments.
!
!     INPUT VARIABLES:
!     IYD - YEAR AND DAY AS YYDDD (day of year from 1 to 365 (or 366))
!     (Year ignored in current model)
!     SEC - UT(SEC)
!     ALT - ALTITUDE(KM)
!     GLAT - GEODETIC LATITUDE(DEG)
!     GLONG - GEODETIC LONGITUDE(DEG)
!     STL - LOCAL APPARENT SOLAR TIME(HRS; see Note below)
!     F107A - 81 day AVERAGE OF F10.7 FLUX (centered on day DDD)
!     F107 - DAILY F10.7 FLUX FOR PREVIOUS DAY
!     AP - MAGNETIC INDEX(DAILY) OR WHEN SW(9)=-1. :
!     - ARRAY CONTAINING:
!     (1) DAILY AP
!     (2) 3 HR AP INDEX FOR CURRENT TIME
!     (3) 3 HR AP INDEX FOR 3 HRS BEFORE CURRENT TIME
!     (4) 3 HR AP INDEX FOR 6 HRS BEFORE CURRENT TIME
!     (5) 3 HR AP INDEX FOR 9 HRS BEFORE CURRENT TIME
!     (6) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 12 TO 33 HRS PRIOR
!     TO CURRENT TIME
!     (7) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 36 TO 57 HRS PRIOR
!     TO CURRENT TIME
!     MASS - MASS NUMBER (ONLY DENSITY FOR SELECTED GAS IS
!     CALCULATED.  MASS 0 IS TEMPERATURE.  MASS 48 FOR ALL.
!     MASS 17 IS Anomalous O ONLY.)
!
!     NOTES ON INPUT VARIABLES:
!     UT, Local Time, and Longitude are used independently in the
!     model and are not of equal importance for every situation.
!     For the most physically realistic calculation these three
!     variables should be consistent (STL=SEC/3600+GLONG/15).
!     The Equation of Time departures from the above formula
!     for apparent local time can be included if available but
!     are of minor importance.
!
!     F107 and F107A values used to generate the model correspond
!     to the 10.7 cm radio flux at the actual distance of the Earth
!     from the Sun rather than the radio flux at 1 AU.
!
!     OUTPUT VARIABLES:
!     D(1) - HE NUMBER DENSITY(CM-3)
!     D(2) - O NUMBER DENSITY(CM-3)
!     D(3) - N2 NUMBER DENSITY(CM-3)
!     D(4) - O2 NUMBER DENSITY(CM-3)
!     D(5) - AR NUMBER DENSITY(CM-3)
!     D(6) - TOTAL MASS DENSITY(GM/CM3) [includes anomalous oxygen]
!     D(7) - H NUMBER DENSITY(CM-3)
!     D(8) - N NUMBER DENSITY(CM-3)
!     D(9) - Anomalous oxygen NUMBER DENSITY(CM-3)
!     T(1) - EXOSPHERIC TEMPERATURE
!     T(2) - TEMPERATURE AT ALT
!
      call gtd7(iyd,sec,alt,glat,glong,stl,f107a,f107,ap,mass,d,t)
!     TOTAL MASS DENSITY
      if ( mass == 48 ) then
        d(6) = 1.66E-24_r8*(d_4*d(1)+d_16*d(2)+28.0_r8*d(3) + &
                            d_32*d(4)+40.0_r8*d(5)+d(7) +     &
                            14.0_r8*d(8)+d_16*d(9))
        if ( imr == 1 ) d(6) = d(6)/1000.0_r8
      end if
    end subroutine gtd7d
!
    subroutine gtd7(iyd,sec,alt,glat,glong,stl,f107a,f107,ap,mass,d,t)
      implicit none
!
      real(dp) , intent(in) :: alt , f107 , f107a , glat , glong , sec , stl
      integer , intent(in) :: iyd , mass
      real(dp) , intent(in) , dimension(7) :: ap
      real(dp) , intent(out) , dimension(9) :: d
      real(dp) , intent(out) , dimension(2) :: t
!
      real(dp) :: altt , dmc , dmr , dm28m , dz28 , tz , v1 , xlat , xmm
      real(dp) , dimension(9) :: ds
      integer :: j , mss
      real(dp) , dimension(2) :: ts
!
      real(dp) , dimension(25) , save :: sv
      integer , save :: mssl
      real(dp) , save :: alast
      real(dp) , parameter :: zmix = 62.5_r8
!
      data alast /99999.0_r8/
      data mssl /-999/
      data sv /25*d_1/
!
!     NRLMSISE-00
!     -----------
!     Neutral Atmosphere Empirical Model from the surface to lower
!     exosphere
!
!     NEW FEATURES:
!     *Extensive satellite drag database used in model generation
!     *Revised O2 (and O) in lower thermosphere
!     *Additional nonlinear solar activity term
!     *"ANOMALOUS OXYGEN" NUMBER DENSITY, OUTPUT D(9)
!     At high altitudes (> 500 km), hot atomic oxygen or ionized
!     oxygen can become appreciable for some ranges of subroutine
!     inputs, thereby affecting drag on satellites and debris. We
!     group these species under the term "anomalous oxygen," since
!     their individual variations are not presently separable with
!     the drag data used to define this model component.
!
!     SUBROUTINES FOR SPECIAL OUTPUTS:
!
!     HIGH ALTITUDE DRAG: EFFECTIVE TOTAL MASS DENSITY
!     (SUBROUTINE GTD7D, OUTPUT D(6))
!     For atmospheric drag calculations at altitudes above 500 km,
!     call SUBROUTINE GTD7D to compute the "effective total mass
!     density" by including contributions from "anomalous oxygen."
!     See "NOTES ON OUTPUT VARIABLES" below on D(6).
!
!     PRESSURE GRID (SUBROUTINE GHP7)
!     See subroutine GHP7 to specify outputs at a pressure level
!     rather than at an altitude.
!
!     OUTPUT IN M-3 and KG/M3:   CALL METERS(.TRUE.)
!
!     INPUT VARIABLES:
!     IYD - YEAR AND DAY AS YYDDD (day of year from 1 to 365 (or 366))
!     (Year ignored in current model)
!     SEC - UT(SEC)
!     ALT - ALTITUDE(KM)
!     GLAT - GEODETIC LATITUDE(DEG)
!     GLONG - GEODETIC LONGITUDE(DEG)
!     STL - LOCAL APPARENT SOLAR TIME(HRS; see Note below)
!     F107A - 81 day AVERAGE OF F10.7 FLUX (centered on day DDD)
!     F107 - DAILY F10.7 FLUX FOR PREVIOUS DAY
!     AP - MAGNETIC INDEX(DAILY) OR WHEN SW(9)=-1. :
!     - ARRAY CONTAINING:
!     (1) DAILY AP
!     (2) 3 HR AP INDEX FOR CURRENT TIME
!     (3) 3 HR AP INDEX FOR 3 HRS BEFORE CURRENT TIME
!     (4) 3 HR AP INDEX FOR 6 HRS BEFORE CURRENT TIME
!     (5) 3 HR AP INDEX FOR 9 HRS BEFORE CURRENT TIME
!     (6) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 12 TO 33 HRS PRIOR
!     TO CURRENT TIME
!     (7) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 36 TO 57 HRS PRIOR
!     TO CURRENT TIME
!     MASS - MASS NUMBER (ONLY DENSITY FOR SELECTED GAS IS
!     CALCULATED.  MASS 0 IS TEMPERATURE.  MASS 48 FOR ALL.
!     MASS 17 IS Anomalous O ONLY.)
!
!     NOTES ON INPUT VARIABLES:
!     UT, Local Time, and Longitude are used independently in the
!     model and are not of equal importance for every situation.
!     For the most physically realistic calculation these three
!     variables should be consistent (STL=SEC/3600+GLONG/15).
!     The Equation of Time departures from the above formula
!     for apparent local time can be included if available but
!     are of minor importance.
!
!     F107 and F107A values used to generate the model correspond
!     to the 10.7 cm radio flux at the actual distance of the Earth
!     from the Sun rather than the radio flux at 1 AU. The following
!     site provides both classes of values:
!     ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/
!
!     F107, F107A, and AP effects are neither large nor well
!     established below 80 km and these parameters should be set to
!     150., 150., and 4. respectively.
!
!     OUTPUT VARIABLES:
!     D(1) - HE NUMBER DENSITY(CM-3)
!     D(2) - O NUMBER DENSITY(CM-3)
!     D(3) - N2 NUMBER DENSITY(CM-3)
!     D(4) - O2 NUMBER DENSITY(CM-3)
!     D(5) - AR NUMBER DENSITY(CM-3)
!     D(6) - TOTAL MASS DENSITY(GM/CM3)
!     D(7) - H NUMBER DENSITY(CM-3)
!     D(8) - N NUMBER DENSITY(CM-3)
!     D(9) - Anomalous oxygen NUMBER DENSITY(CM-3)
!     T(1) - EXOSPHERIC TEMPERATURE
!     T(2) - TEMPERATURE AT ALT
!
!     NOTES ON OUTPUT VARIABLES:
!     TO GET OUTPUT IN M-3 and KG/M3:   CALL METERS(.TRUE.)
!
!     O, H, and N are set to zero below 72.5 km
!
!     T(1), Exospheric temperature, is set to global average for
!     altitudes below 120 km. The 120 km gradient is left at global
!     average value for altitudes below 72 km.
!
!     D(6), TOTAL MASS DENSITY, is NOT the same for subroutines GTD7
!     and GTD7D
!
!     SUBROUTINE GTD7 -- D(6) is the sum of the mass densities of the
!     species labeled by indices 1-5 and 7-8 in output variable D.
!     This includes He, O, N2, O2, Ar, H, and N but does NOT include
!     anomalous oxygen (species index 9).
!
!     SUBROUTINE GTD7D -- D(6) is the "effective total mass density
!     for drag" and is the sum of the mass densities of all species
!     in this model, INCLUDING anomalous oxygen.
!
!     SWITCHES: The following is for test and special purposes:
!
!     TO TURN ON AND OFF PARTICULAR VARIATIONS CALL TSELEC(SW),
!     WHERE SW IS A 25 ELEMENT ARRAY CONTAINING 0. FOR OFF, 1.
!     FOR ON, OR 2. FOR MAIN EFFECTS OFF BUT CROSS TERMS ON
!     FOR THE FOLLOWING VARIATIONS
!     1 - F10.7 EFFECT ON MEAN  2 - TIME INDEPENDENT
!     3 - SYMMETRICAL ANNUAL    4 - SYMMETRICAL SEMIANNUAL
!     5 - ASYMMETRICAL ANNUAL   6 - ASYMMETRICAL SEMIANNUAL
!     7 - DIURNAL               8 - SEMIDIURNAL
!     9 - DAILY AP             10 - ALL UT/LONG EFFECTS
!     11 - LONGITUDINAL         12 - UT AND MIXED UT/LONG
!     13 - MIXED AP/UT/LONG     14 - TERDIURNAL
!     15 - DEPARTURES FROM DIFFUSIVE EQUILIBRIUM
!     16 - ALL TINF VAR         17 - ALL TLB VAR
!     18 - ALL TN1 VAR           19 - ALL S VAR
!     20 - ALL TN2 VAR           21 - ALL NLB VAR
!     22 - ALL TN3 VAR           23 - TURBO SCALE HEIGHT VAR
!
!     To get current values of SW: CALL TRETRV(SW)
!
      if ( isw /= 64999 ) call tselec(sv)
!
      ! Test for changed input
      v1 = vtst7(iyd,sec,glat,glong,stl,f107a,f107,ap,1)
      !
      ! Latitude variation of gravity (none for SW(2)=0)
      xlat = glat
      if ( dabs(sw(2)) < nearzero ) xlat = d_045
      call glatf(xlat)
!
      xmm = pdm(5,3)
!
      ! THERMOSPHERE/MESOSPHERE (above ZN2(1))
      altt = dmax1(alt,zn2(1))
      mss = mass
      ! Only calculate N2 in thermosphere if alt in mixed region
      if ( alt < zmix .and. mass > 0 ) mss = 28
      ! Only calculate thermosphere if input parameters changed
      ! or altitude above ZN2(1) in mesosphere
      if ( dabs(v1-d_1) < nearzero .or. alt > zn2(1) .or. &
           alast > zn2(1) .or. mss /= mssl ) then
        call gts7(iyd,sec,altt,glat,glong,stl,f107a,f107,ap,mss,ds,ts)
        dm28m = dm28
!       metric adjustment
        if ( imr == 1 ) dm28m = dm28*1.E6_r8
        mssl = mss
      end if
      t(1) = ts(1)
      t(2) = ts(2)
      if ( alt >= zn2(1) ) then
        do j = 1 , 9
          d(j) = ds(j)
        end do
        return
      end if
      !
      ! LOWER MESOSPHERE/UPPER STRATOSPHERE [between ZN3(1) and ZN2(1)]
      ! Temperature at nodes and gradients at end nodes
      ! Inverse temperature a linear function of spherical harmonics
      ! Only calculate nodes if input changed
      !
      if ( dabs(v1-d_1) < nearzero .or. alast >= zn2(1) ) then
        tgn2(1) = tgn1(2)
        tn2(1) = tn1(5)
        tn2(2) = pma(1,1)*pavgm(1)/(d_1-sw(20)*glob7s(glong,pma(:,1)))
        tn2(3) = pma(1,2)*pavgm(2)/(d_1-sw(20)*glob7s(glong,pma(:,2)))
        tn2(4) = pma(1,3)*pavgm(3)/(d_1-sw(20)*sw(22)*glob7s(glong,pma(:,3)))
        tgn2(2) = pavgm(9)*pma(1,10)                                    &
                 *(d_1+sw(20)*sw(22)*glob7s(glong,pma(:,10)))*tn2(4)*tn2(4)   &
                 /(pma(1,3)*pavgm(3))**d_2
        tn3(1) = tn2(4)
      end if
      !
      ! LOWER STRATOSPHERE AND TROPOSPHERE [below ZN3(1)]
      ! Temperature at nodes and gradients at end nodes
      ! Inverse temperature a linear function of spherical harmonics
      ! Only calculate nodes if input changed
      !
      if ( alt < zn3(1) ) then
        if ( dabs(v1-d_1) < nearzero .or. alast >= zn3(1) ) then
          tgn3(1) = tgn2(2)
          tn3(2) = pma(1,4)*pavgm(4)/(d_1-sw(22)*glob7s(glong,pma(:,4)))
          tn3(3) = pma(1,5)*pavgm(5)/(d_1-sw(22)*glob7s(glong,pma(:,5)))
          tn3(4) = pma(1,6)*pavgm(6)/(d_1-sw(22)*glob7s(glong,pma(:,6)))
          tn3(5) = pma(1,7)*pavgm(7)/(d_1-sw(22)*glob7s(glong,pma(:,7)))
          tgn3(2) = pma(1,8)*pavgm(8)*(d_1+sw(22)*glob7s(glong,pma(:,8)))  &
                   *tn3(5)*tn3(5)/(pma(1,7)*pavgm(7))**d_2
        end if
      end if

      if ( mass == 0 ) then
        dd = densm(alt,d_1,d_0,tz)
        t(2) = tz
      else
        ! LINEAR TRANSITION TO FULL MIXING BELOW ZN2(1)
        dmc = d_0
        if ( alt > zmix ) dmc = d_1 - (zn2(1)-alt)/(zn2(1)-zmix)
        dz28 = ds(3)
        ! N2 DENSITY
        dmr = ds(3)/dm28m - d_1
        d(3) = densm(alt,dm28m,xmm,tz)
        d(3) = d(3)*(d_1+dmr*dmc)
        ! HE DENSITY
        d(1) = d_0
        if ( mass == 4 .or. mass == 48 ) then
          dmr = ds(1)/(dz28*pdm(2,1)) - d_1
          d(1) = d(3)*pdm(2,1)*(d_1+dmr*dmc)
        end if
        ! O DENSITY
        d(2) = d_0
        d(9) = d_0
        ! O2 DENSITY
        d(4) = d_0
        if ( mass == 32 .or. mass == 48 ) then
          dmr = ds(4)/(dz28*pdm(2,4)) - d_1
          d(4) = d(3)*pdm(2,4)*(d_1+dmr*dmc)
        end if
        ! AR DENSITY
        d(5) = d_0
        if ( mass == 40 .or. mass == 48 ) then
          dmr = ds(5)/(dz28*pdm(2,5)) - d_1
          d(5) = d(3)*pdm(2,5)*(d_1+dmr*dmc)
        end if
        ! HYDROGEN DENSITY
        d(7) = d_0
        ! ATOMIC NITROGEN DENSITY
        d(8) = d_0
        !
        ! TOTAL MASS DENSITY
        !
        if ( mass == 48 ) then
          d(6) = 1.66E-24_r8*(d_4*d(1)+d_16*d(2)+28.0_r8*d(3) + &
                              d_32*d(4)+40.0_r8*d(5)+d(7)+14.0_r8*d(8))
          if ( imr == 1 ) d(6) = d(6)/1000.0_r8
        end if
        t(2) = tz
      end if
      alast = alt
    end subroutine gtd7
!
    subroutine gts7(iyd,sec,alt,glat,glong,stl,f107a,f107,ap,mass,d,t)
      implicit none
!
      real(dp) , intent(in) :: alt , f107 , f107a , glat , glong , sec , stl
      integer , intent(in) :: iyd , mass
      real(dp) , intent(in) , dimension(:) :: ap
      real(dp) , intent(inout) , dimension(9) :: d
      real(dp) , intent(inout) , dimension(2) :: t
!
      real(dp) :: b01 , b04 , b14 , b16 , b28 , b32 , b40 ,     &
              day , db16h , ddum , g1 , g14 , g16 , g16h , g28 ,        &
              g32 , g4 , g40 , hc01 , hc04 , hc14 , hc16 , hc216 ,      &
              hc32 , hc40 , hcc01 , hcc14 , hcc16 , hcc232 , hcc32 ,    &
              rc01 , rc14 , rc16 , rc32 , t2 , tho , tz , v2 ,   &
              xmd , xmm , yrd , z , zc01 , zc04 , zc14 , zc16 , zc32 ,  &
              zc40 , zcc01 , zcc14 , zcc16 , zcc32 , zh01 , zh04 ,      &
              zh14 , zh16 , zh28 , zh32 , zh40 , zhf , zhm01 , zhm04 ,  &
              zhm14 , zhm16 , zhm28 , zhm32 , zhm40 , zmho , zsho , zsht
      integer :: i , j
      real(dp) , save :: alast
!
!     Thermospheric portion of NRLMSISE-00
!     See GTD7 for more extensive comments
!
!     OUTPUT IN M-3 and KG/M3:   CALL METERS(.TRUE.)
!
!     INPUT VARIABLES:
!     IYD - YEAR AND DAY AS YYDDD (day of year from 1 to 365 (or 366))
!     (Year ignored in current model)
!     SEC - UT(SEC)
!     ALT - ALTITUDE(KM) (>72.5 km)
!     GLAT - GEODETIC LATITUDE(DEG)
!     GLONG - GEODETIC LONGITUDE(DEG)
!     STL - LOCAL APPARENT SOLAR TIME(HRS; see Note below)
!     F107A - 81 day AVERAGE OF F10.7 FLUX (centered on day DDD)
!     F107 - DAILY F10.7 FLUX FOR PREVIOUS DAY
!     AP - MAGNETIC INDEX(DAILY) OR WHEN SW(9)=-1. :
!     - ARRAY CONTAINING:
!     (1) DAILY AP
!     (2) 3 HR AP INDEX FOR CURRENT TIME
!     (3) 3 HR AP INDEX FOR 3 HRS BEFORE CURRENT TIME
!     (4) 3 HR AP INDEX FOR 6 HRS BEFORE CURRENT TIME
!     (5) 3 HR AP INDEX FOR 9 HRS BEFORE CURRENT TIME
!     (6) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 12 TO 33 HRS PRIOR
!     TO CURRENT TIME
!     (7) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 36 TO 57 HRS PRIOR
!     TO CURRENT TIME
!     MASS - MASS NUMBER (ONLY DENSITY FOR SELECTED GAS IS
!     CALCULATED.  MASS 0 IS TEMPERATURE.  MASS 48 FOR ALL.
!     MASS 17 IS Anomalous O ONLY.)
!
!     NOTES ON INPUT VARIABLES:
!     UT, Local Time, and Longitude are used independently in the
!     model and are not of equal importance for every situation.
!     For the most physically realistic calculation these three
!     variables should be consistent (STL=SEC/3600+GLONG/15).
!     The Equation of Time departures from the above formula
!     for apparent local time can be included if available but
!     are of minor importance.
!
!     F107 and F107A values used to generate the model correspond
!     to the 10.7 cm radio flux at the actual distance of the Earth
!     from the Sun rather than the radio flux at 1 AU. The following
!     site provides both classes of values:
!     ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/
!
!     F107, F107A, and AP effects are neither large nor well
!     established below 80 km and these parameters should be set to
!     150., 150., and 4. respectively.
!
!     OUTPUT VARIABLES:
!     D(1) - HE NUMBER DENSITY(CM-3)
!     D(2) - O NUMBER DENSITY(CM-3)
!     D(3) - N2 NUMBER DENSITY(CM-3)
!     D(4) - O2 NUMBER DENSITY(CM-3)
!     D(5) - AR NUMBER DENSITY(CM-3)
!     D(6) - TOTAL MASS DENSITY(GM/CM3) [Anomalous O NOT included]
!     D(7) - H NUMBER DENSITY(CM-3)
!     D(8) - N NUMBER DENSITY(CM-3)
!     D(9) - Anomalous oxygen NUMBER DENSITY(CM-3)
!     T(1) - EXOSPHERIC TEMPERATURE
!     T(2) - TEMPERATURE AT ALT
!
      data alast / -999.0_r8/
!
      ! Test for changed input
      v2 = vtst7(iyd,sec,glat,glong,stl,f107a,f107,ap,2)
!
      yrd = iyd
      za = pdl(16,2)
      zn1(1) = za
      do j = 1 , 9
        d(j) = d_0
      end do
      ! TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1)
      if ( alt > zn1(1) ) then
        if ( dabs(v2-d_1) < nearzero .or. alast <= zn1(1) ) then
          tinf = ptm(1)*pt(1)*(d_1+sw(16)* &
             globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pt))
        end if
      else
        tinf = ptm(1)*pt(1)
      end if
      t(1) = tinf
      ! GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5)
      if ( alt > zn1(5) ) then
        if ( dabs(v2-d_1) < nearzero .or. alast <= zn1(5) ) then
          xg0 = ptm(4)*ps(1)*(d_1+sw(19) * &
             globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,ps))
        end if
      else
        xg0 = ptm(4)*ps(1)
      end if
      ! Calculate these temperatures only if input changed
      if ( dabs(v2-d_1) < nearzero .or. alt < 300.0_r8 ) then
        tlb = ptm(2)*(d_1+sw(17)* &
          globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(:,4)))*pd(1,4)
      end if
      s = xg0/(tinf-tlb)
      ! Lower thermosphere temp variations not significant for
      ! density above 300 km
      if ( alt >= 300.0_r8 ) then
        tn1(2) = ptm(7)*ptl(1,1)
        tn1(3) = ptm(3)*ptl(1,2)
        tn1(4) = ptm(8)*ptl(1,3)
        tn1(5) = ptm(5)*ptl(1,4)
        tgn1(2) = ptm(9)*pma(1,9)*tn1(5)*tn1(5)/(ptm(5)*ptl(1,4))**d_2
      else if ( dabs(v2-d_1) < nearzero .or. alast >= 300.0_r8 ) then
        tn1(2) = ptm(7)*ptl(1,1)/(d_1-sw(18)*glob7s(glong,ptl(:,1)))
        tn1(3) = ptm(3)*ptl(1,2)/(d_1-sw(18)*glob7s(glong,ptl(:,2)))
        tn1(4) = ptm(8)*ptl(1,3)/(d_1-sw(18)*glob7s(glong,ptl(:,3)))
        tn1(5) = ptm(5)*ptl(1,4)/(d_1-sw(18)*sw(20)*glob7s(glong,ptl(:,4)))
        tgn1(2) = ptm(9)*pma(1,9)*(d_1+sw(18)*sw(20)*glob7s(glong,pma(:,9)))   &
                & *tn1(5)*tn1(5)/(ptm(5)*ptl(1,4))**d_2
      end if
!
      z0 = zn1(4)
      t0 = tn1(4)
      tr12 = d_1
!
      if ( mass == 0 ) go to 700
!
      ! N2 variation factor at Zlb
      g28 = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(:,3))
      day = dmod(yrd,1000.0_r8)
      ! VARIATION OF TURBOPAUSE HEIGHT
      zhf = pdl(25,2)*(d_1+sw(5) * &
            pdl(25,1)*dsin(deg2rad*glat)*dcos(dr*(day-pt(14))))
      yrd = iyd
      t(1) = tinf
      xmm = pdm(5,3)
      z = alt
!
      do j = 1 , 11
        if ( mass == mt(j) ) go to 100
      end do
!
      go to 800
!
 100  continue
!
      if ( z <= altl(6) .or. mass == 28 .or. mass == 48 ) then
        !
        !       **** N2 DENSITY ****
        !
        ! Diffusive density at Zlb
        db28 = pdm(1,3)*dexp(g28)*pd(1,3)
        ! Diffusive density at Alt
        d(3) = densu(z,db28,tinf,tlb,28.0_r8,alpha(3),t(2),ptm(6),s)
        dd = d(3)
        ! Turbopause
        zh28 = pdm(3,3)*zhf
        zhm28 = pdm(4,3)*pdl(6,2)
        xmd = 28.0_r8 - xmm
        ! Mixed density at Zlb
        b28 = densu(zh28,db28,tinf,tlb,xmd,alpha(3)-d_1,tz,ptm(6),s)
        if ( z <= altl(3) .and. dabs(sw(15)) > d_0 ) then
          ! Mixed density at Alt
          dm28 = densu(z,b28,tinf,tlb,xmm,alpha(3),tz,ptm(6),s)
          ! Net density at Alt
          d(3) = dnet(d(3),dm28,zhm28,xmm,28.0_r8)
        end if
      end if
      if ( j == 2 ) go to 700
      if ( j == 4 .or. j == 9 ) then
      else if ( j == 5 ) then
        go to 800
      else if ( j == 6 ) then
        go to 200
      else if ( j == 7 ) then
        go to 300
      else if ( j == 8 ) then
        go to 400
      else if ( j == 10 ) then
        go to 500
      else if ( j == 11 ) then
        go to 600
      else
        !
        !       **** HE DENSITY ****
        !
        ! Density variation factor at Zlb
        g4 = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(:,1))
        ! Diffusive density at Zlb
        db04 = pdm(1,1)*dexp(g4)*pd(1,1)
        ! Diffusive density at Alt
        d(1) = densu(z,db04,tinf,tlb,d_4,alpha(1),t(2),ptm(6),s)
        dd = d(1)
        if ( z <= altl(1) .and. dabs(sw(15)) > d_0 ) then
          ! Turbopause
          zh04 = pdm(3,1)
          ! Mixed density at Zlb
          b04 = densu(zh04,db04,tinf,tlb,d_4-xmm,alpha(1)-d_1,t(2),ptm(6),s)
          ! Mixed density at Alt
          dm04 = densu(z,b04,tinf,tlb,xmm,d_0,t(2),ptm(6),s)
          zhm04 = zhm28
          ! Net density at Alt
          d(1) = dnet(d(1),dm04,zhm04,xmm,d_4)
          ! Correction to specified mixing ratio at ground
          rl = dlog(b28*pdm(2,1)/b04)
          zc04 = pdm(5,1)*pdl(1,2)
          hc04 = pdm(6,1)*pdl(2,2)
          ! Net density corrected at Alt
          d(1) = d(1)*ccor(z,rl,hc04,zc04)
        end if
        if ( mass /= 48 ) go to 800
      end if
      !
      !     **** O DENSITY ****
      !
      ! Density variation factor at Zlb
      g16 = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(:,2))
      ! Diffusive density at Zlb
      db16 = pdm(1,2)*dexp(g16)*pd(1,2)
      ! Diffusive density at Alt
      d(2) = densu(z,db16,tinf,tlb,d_16,alpha(2),t(2),ptm(6),s)
      dd = d(2)
      if ( z <= altl(2) .and. dabs(sw(15)) > d_0 ) then
        ! Corrected from PDM(3,1) to PDM(3,2)  12/2/85
        ! Turbopause
        zh16 = pdm(3,2)
        ! Mixed density at Zlb
        b16 = densu(zh16,db16,tinf,tlb,d_16-xmm,alpha(2)-d_1,t(2),ptm(6),s)
        ! Mixed density at Alt
        dm16 = densu(z,b16,tinf,tlb,xmm,d_0,t(2),ptm(6),s)
        zhm16 = zhm28
        ! Net density at Alt
        d(2) = dnet(d(2),dm16,zhm16,xmm,d_16)
        ! 3/16/99 Change form to match O2 departure from diff equil near
        ! 150 km and add dependence on F10.7
        ! RL=ALOG(B28*PDM(2,2)*ABS(PDL(17,2))/B16)
        rl = pdm(2,2)*pdl(17,2)*(d_1+sw(1)*pdl(24,1)*(f107a-150.0_r8))
        hc16 = pdm(6,2)*pdl(4,2)
        zc16 = pdm(5,2)*pdl(3,2)
        hc216 = pdm(6,2)*pdl(5,2)
        d(2) = d(2)*ccor2(z,rl,hc16,zc16,hc216)
        ! Chemistry correction
        hcc16 = pdm(8,2)*pdl(14,2)
        zcc16 = pdm(7,2)*pdl(13,2)
        rc16 = pdm(4,2)*pdl(15,2)
        ! Net density corrected at Alt
        d(2) = d(2)*ccor(z,rc16,hcc16,zcc16)
      end if
      if ( mass /= 48 .and. mass /= 49 ) go to 800
      !
      !     **** O2 DENSITY ****
      !
      ! Density variation factor at Zlb
 200  continue
      g32 = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(:,5))
      ! Diffusive density at Zlb
      db32 = pdm(1,4)*dexp(g32)*pd(1,5)
      ! Diffusive density at Alt
      d(4) = densu(z,db32,tinf,tlb,d_32,alpha(4),t(2),ptm(6),s)
      if ( mass == 49 ) then
        dd = dd + d_2*d(4)
      else
        dd = d(4)
      end if
      if ( dabs(sw(15)) > d_0 ) then
        if ( z <= altl(4) ) then
          ! Turbopause
          zh32 = pdm(3,4)
          ! Mixed density at Zlb
          b32 = densu(zh32,db32,tinf,tlb,d_32-xmm,alpha(4)-d_1, t(2),ptm(6),s)
          ! Mixed density at Alt
          dm32 = densu(z,b32,tinf,tlb,xmm,d_0,t(2),ptm(6),s)
          zhm32 = zhm28
          ! Net density at Alt
          d(4) = dnet(d(4),dm32,zhm32,xmm,d_32)
          ! Correction to specified mixing ratio at ground
          rl = dlog(b28*pdm(2,4)/b32)
          hc32 = pdm(6,4)*pdl(8,2)
          zc32 = pdm(5,4)*pdl(7,2)
          d(4) = d(4)*ccor(z,rl,hc32,zc32)
        end if
        ! Correction for general departure from diffusive equilibrium
        ! above Zlb
        hcc32 = pdm(8,4)*pdl(23,2)
        hcc232 = pdm(8,4)*pdl(23,1)
        zcc32 = pdm(7,4)*pdl(22,2)
        rc32 = pdm(4,4)*pdl(24,2)*(d_1+sw(1)*pdl(24,1)*(f107a-150.0_r8))
        ! Net density corrected at Alt
        d(4) = d(4)*ccor2(z,rc32,hcc32,zcc32,hcc232)
      end if
      if ( mass /= 48 ) go to 800
      !
      !     **** AR DENSITY ****
      !
      ! Density variation factor at Zlb
 300  continue
      g40 = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(:,6))
      ! Diffusive density at Zlb
      db40 = pdm(1,5)*dexp(g40)*pd(1,6)
      ! Diffusive density at Alt
      d(5) = densu(z,db40,tinf,tlb,40.0_r8,alpha(5),t(2),ptm(6),s)
      dd = d(5)
      if ( z <= altl(5) .and. dabs(sw(15)) > d_0 ) then
        ! Turbopause
        zh40 = pdm(3,5)
        ! Mixed density at Zlb
        b40 = densu(zh40,db40,tinf,tlb,40.0_r8-xmm,alpha(5)-d_1,t(2),ptm(6),s)
        ! Mixed density at Alt
        dm40 = densu(z,b40,tinf,tlb,xmm,d_0,t(2),ptm(6),s)
        zhm40 = zhm28
        ! Net density at Alt
        d(5) = dnet(d(5),dm40,zhm40,xmm,40.0_r8)
        ! Correction to specified mixing ratio at ground
        rl = dlog(b28*pdm(2,5)/b40)
        hc40 = pdm(6,5)*pdl(10,2)
        zc40 = pdm(5,5)*pdl(9,2)
        ! Net density corrected at Alt
        d(5) = d(5)*ccor(z,rl,hc40,zc40)
      end if
      if ( mass /= 48 ) go to 800
      !
      !     **** HYDROGEN DENSITY ****
      !
      ! Density variation factor at Zlb
 400  continue
      g1 = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(:,7))
      ! Diffusive density at Zlb
      db01 = pdm(1,6)*dexp(g1)*pd(1,7)
      ! Diffusive density at Alt
      d(7) = densu(z,db01,tinf,tlb,d_1,alpha(7),t(2),ptm(6),s)
      dd = d(7)
      if ( z <= altl(7) .and. dabs(sw(15)) > d_0 ) then
        ! Turbopause
        zh01 = pdm(3,6)
        ! Mixed density at Zlb
        b01 = densu(zh01,db01,tinf,tlb,d_1-xmm,alpha(7)-d_1,t(2),ptm(6),s)
        ! Mixed density at Alt
        dm01 = densu(z,b01,tinf,tlb,xmm,d_0,t(2),ptm(6),s)
        zhm01 = zhm28
        ! Net density at Alt
        d(7) = dnet(d(7),dm01,zhm01,xmm,d_1)
        ! Correction to specified mixing ratio at ground
        rl = dlog(b28*pdm(2,6)*dabs(pdl(18,2))/b01)
        hc01 = pdm(6,6)*pdl(12,2)
        zc01 = pdm(5,6)*pdl(11,2)
        d(7) = d(7)*ccor(z,rl,hc01,zc01)
        ! Chemistry correction
        hcc01 = pdm(8,6)*pdl(20,2)
        zcc01 = pdm(7,6)*pdl(19,2)
        rc01 = pdm(4,6)*pdl(21,2)
        ! Net density corrected at Alt
        d(7) = d(7)*ccor(z,rc01,hcc01,zcc01)
      end if
      if ( mass /= 48 ) go to 800
      !
      !     **** ATOMIC NITROGEN DENSITY ****
      !
      ! Density variation factor at Zlb
 500  continue
      g14 = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(:,8))
      ! Diffusive density at Zlb
      db14 = pdm(1,7)*dexp(g14)*pd(1,8)
      ! Diffusive density at Alt
      d(8) = densu(z,db14,tinf,tlb,14.0_r8,alpha(8),t(2),ptm(6),s)
      dd = d(8)
      if ( z <= altl(8) .and. dabs(sw(15)) > d_0 ) then
        ! Turbopause
        zh14 = pdm(3,7)
        ! Mixed density at Zlb
        b14 = densu(zh14,db14,tinf,tlb,14.0_r8-xmm,alpha(8)-d_1,t(2),ptm(6),s)
        ! Mixed density at Alt
        dm14 = densu(z,b14,tinf,tlb,xmm,d_0,t(2),ptm(6),s)
        zhm14 = zhm28
        ! Net density at Alt
        d(8) = dnet(d(8),dm14,zhm14,xmm,14.0_r8)
        ! Correction to specified mixing ratio at ground
        rl = dlog(b28*pdm(2,7)*dabs(pdl(3,1))/b14)
        hc14 = pdm(6,7)*pdl(2,1)
        zc14 = pdm(5,7)*pdl(1,1)
        d(8) = d(8)*ccor(z,rl,hc14,zc14)
        ! Chemistry correction
        hcc14 = pdm(8,7)*pdl(5,1)
        zcc14 = pdm(7,7)*pdl(4,1)
        rc14 = pdm(4,7)*pdl(6,1)
        ! Net density corrected at Alt
        d(8) = d(8)*ccor(z,rc14,hcc14,zcc14)
      end if
      if ( mass /= 48 ) go to 800
      !
      !     **** Anomalous OXYGEN DENSITY ****
      !
 600  continue
      g16h = sw(21)*globe7(yrd,sec,glat,glong,stl,f107a,f107,ap,pd(:,9))
      db16h = pdm(1,8)*dexp(g16h)*pd(1,9)
      tho = pdm(10,8)*pdl(7,1)
      dd = densu(z,db16h,tho,tho,d_16,alpha(9),t2,ptm(6),s)
      zsht = pdm(6,8)
      zmho = pdm(5,8)
      zsho = scalh(zmho,d_16,tho)
      d(9) = dd*dexp(-zsht/zsho*(dexp(-(z-zmho)/zsht)-d_1))
      if ( mass == 48 ) then
        !
        ! TOTAL MASS DENSITY
        !
        d(6) = 1.66E-24_r8*(d_4*d(1)+d_16*d(2)+28.0_r8*d(3) + &
                            d_32*d(4)+40.0_r8*d(5)+d(7)+14.0_r8*d(8))
        db48 = 1.66E-24*(d_4*db04+d_16*db16+28.0_r8*db28+d_32*db32 + &
                         40.0_r8*db40+db01+14.0_r8*db14)
      end if
      go to 800
      ! TEMPERATURE AT ALTITUDE
 700  continue
      z = dabs(alt)
      ddum = densu(z,d_1,tinf,tlb,d_0,d_0,t(2),ptm(6),s)
      ! ADJUST DENSITIES FROM CGS TO KGM
 800  continue
      if ( imr == 1 ) then
        do i = 1 , 9
          d(i) = d(i)*1.E6_r8
        end do
        d(6) = d(6)/1000.0_r8
      end if
      alast = alt
      return

      contains
      !
      !   Calculate scale height (km)
      !
      real(dp) function scalh(alt,xm,temp)
        implicit none
        real(dp) , intent(in) :: alt , temp , xm
        real(dp) :: g
        g = gsurf/(d_1+alt/re)**d_2
        scalh = r100gas*temp/(g*xm)
      end function scalh
    end subroutine gts7
!
    real(dp) function vtst7(iyd,sec,glat,glong,stl,f107a,f107,ap,ic)
      implicit none
!
      real(dp) , intent(in) :: f107 , f107a , glat , glong , sec , stl
      integer , intent(in) :: ic , iyd
      real(dp) , intent(in) , dimension(7) :: ap
!
      real(dp) , dimension(7,2) , save :: apl
      real(dp) , dimension(2) , save :: fal , fl , glatl , gll , secl , stll
      integer , dimension(2) , save :: iydl
      real(dp) , dimension(25,2) , save :: swcl , swl
      integer :: i
      logical :: ldc
!
!     Test if geophysical variables or switches changed and save
!     Return 0 if unchanged and 1 if changed
!
      data iydl  /2* -999/
      data secl  /2* -999.0_r8/
      data glatl /2* -999.0_r8/
      data gll   /2* -999.0_r8/
      data stll  /2* -999.0_r8/
      data fal   /2* -999.0_r8/
      data fl    /2* -999.0_r8/
      data apl   /14* -999.0_r8/
      data swl   /50* -999.0_r8/
      data swcl  /50* -999.0_r8/
      vtst7 = d_0
      ldc = .false.

      if ( iyd /= iydl(ic) ) then
        ldc = .true.
      else
        if ( dabs(sec-secl(ic)) > nearzero ) then
          ldc = .true.
        else
          if ( dabs(glat-glatl(ic)) > nearzero ) then
            ldc = .true.
          else
            if ( dabs(glong-gll(ic)) > nearzero ) then
              ldc = .true.
            else
              if ( dabs(stl-stll(ic)) > nearzero ) then
                ldc = .true.
              else
                if ( dabs(f107a-fal(ic)) > nearzero ) then
                  ldc = .true.
                else
                  if ( dabs(f107-fl(ic)) > nearzero ) then
                    ldc = .true.
                  else
                    do i = 1 , 7
                      if ( dabs(ap(i)-apl(i,ic)) > nearzero ) then
                        ldc = .true.
                        exit
                      end if
                    end do
                    do i = 1 , 25
                      if ( dabs(sw(i)-swl(i,ic)) > nearzero ) then
                        ldc = .true.
                        exit
                      end if
                      if ( dabs(swc(i)-swcl(i,ic)) > nearzero ) then
                        ldc = .true.
                        exit
                      end if
                    end do
                  end if
                end if
              end if
            end if
          end if
        end if
      end if

      if ( ldc ) then
        vtst7 = d_1
        iydl(ic) = iyd
        secl(ic) = sec
        glatl(ic) = glat
        gll(ic) = glong
        stll(ic) = stl
        fal(ic) = f107a
        fl(ic) = f107
        do i = 1 , 7
          apl(i,ic) = ap(i)
        end do
        do i = 1 , 25
          swl(i,ic) = sw(i)
          swcl(i,ic) = swc(i)
        end do
      end if

    end function vtst7
!
end module physics_msis
